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Abstract 


Event Horizon Telescope (EHT) observations at 230 GHz have now imaged polarized emission around the 
supermassive black hole in M87 on event-horizon scales. This polarized synchrotron radiation probes the structure 
of magnetic fields and the plasma properties near the black hole. Here we compare the resolved polarization 
structure observed by the EHT, along with simultaneous unresolved observations with the Atacama Large 
Millimeter/submillimeter Array, to expectations from theoretical models. The low fractional linear polarization in 
the resolved image suggests that the polarization is scrambled on scales smaller than the EHT beam, which we 
attribute to Faraday rotation internal to the emission region. We estimate the average density л, 10^7 cm ^, 
magnetic field strength B ~ 1—30 С, and electron temperature T, ~ (1-12) х 10'? K of the radiating plasma in a 
simple one-zone emission model. We show that the net azimuthal linear polarization pattern may result from 
organized, poloidal magnetic fields in the emission region. In a quantitative comparison with a large library of 
simulated polarimetric images from general relativistic magnetohydrodynamic (GRMHD) simulations, we identify 
a subset of physical models that can explain critical features of the polarimetric EHT observations while producing 
a relativistic jet of sufficient power. The consistent GRMHD models are all of magnetically arrested accretion 
disks, where near-horizon magnetic fields are dynamically important. We use the models to infer a mass accretion 
rate onto the black hole in M87 of (3-20) x 10 * Mo yr. 


Unified Astronomy Thesaurus concepts: Accretion (14); Black holes (162); Event horizons (479); Jets (870); Kerr 
black holes (886); Magnetic fields (994); Magnetohydrodynamics (1964); Plasma astrophysics (1261); Polarimetry 


CrossMark 


(1278); Radiative transfer (1335); Radio jets (1347); Relativistic jets (1390) 


1. Introduction 


The Event Horizon Telescope (EHT) Collaboration has 
recently published total intensity images of event-horizon-scale 
emission around the supermassive black hole in the core of the 
M87 galaxy (M87*; Event Horizon Telescope Collaboration et al. 
2019a, 2019b, 2019c, 2019d, hereafter EHTC I, EHTC II, 
EHTC III, EHTC IV). The data reveal a 42 + 3 pas diameter ring- 
like structure that is broadly consistent with the shadow of a black 
hole as predicted by Einstein's Theory of General Relativity 
(Event Horizon Telescope Collaboration et al. 2019e, 2019f; 
hereafter EHTC V, EHTC VI). The brightness temperature of the 
ring at 230 GHz (10/9 K) is naturally explained by synchrotron 
emission from relativistic electrons gyrating around magnetic field 
lines. The ring brightness asymmetry results from light bending 
and Doppler beaming due to relativistic rotation of the matter 
around the black hole. 

М87* is best known for launching a kpc-scale FR-I type 
relativistic jet, whose kinetic power is estimated to be 
104274 erg s^! (e.g., Stawarz et al. 2006; de Gasperin et al. 
2012). The structure of the relativistic jet has been resolved and 


126 NASA Hubble Fellowship Program, Einstein Fellow. 


127 BACOA Fellow. 
128 UKRI Stephen Hawking Fellow. 


Original content from this work may be used under the terms 

BY of the Creative Commons Attribution 4.0 licence. Any further 
distribution of this work must maintain attribution to the author(s) and the title 
of the work, journal citation and DOI. 


studied at radio to X-ray wavelengths (e.g., Di Matteo et al. 
2003; Harris et al. 2009; Kim et al. 2018; Walker et al. 2018). 

The published ЕНТ image of M87* together with multi- 
wavelength observations are consistent with the picture that the 
supermassive black hole in M87 is surrounded by a 
relativistically hot, magnetized plasma (Rees et al. 1982; 
Narayan & Yi 1995; Narayan et al. 1995; Yuan & 
Narayan 2014; Reynolds et al. 1996; Yuan et al. 2002; Di 
Matteo et al. 2003). However, it is not clear whether the 
compact ring emission is produced by plasma that is inflowing 
(in a thick accretion flow), outflowing (at the jet base or in a 
wind), or both. Furthermore, the total intensity EHT observa- 
tions also could not constrain the structure of magnetic fields in 
the observed emission region. In order to find out which 
physical scenario is realized in M87", additional information is 
necessary. 

Event Horizon Telescope Collaboration et al. (2021, 
hereafter EHTC VIT) reports new results from the polarimetric 
ЕНТ 2017 observations of M87". The polarimetric images of 
M87" are reproduced in Figure 1. These images reveal that a 
significant fraction of the ring emission is linearly polarized, as 
expected for synchrotron radiation. The EHT polarimetric 
measurements are consistent with unresolved observations of 
the radio core at the same frequency with the Submillimeter 
Array (SMA; Kuo et al. 2014) and the Atacama Large 
Millimeter/submillimeter Array (ALMA; Goddi et al. 2021). 
They also provide a detailed view of the polarized emission 
region on event-horizon scales near the black hole. Polarized 
synchrotron radiation traces the underlying magnetic field 
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Figure 1. Top panel: 2017 April 11 fiducial polarimetric image of 
M87* from EHTC VII. The gray scale encodes the total intensity, and ticks 
illustrate the degree and direction of linear polarization. The tick color indicates 
the amplitude of the fractional linear polarization, the tick length is proportional 


to |P] = JQ? + U7’, and the tick direction indicates the electric-vector 
position angle (EVPA). Polarization ticks are displayed only in regions where 
T > 10% Zma and |P| > 20%|Plmax. Bottom row: polarimetric images of 
M87* taken on different days. 


configuration and magnetized plasma properties along the line 
of sight (Bromley et al. 2001; Broderick & Loeb 2009; 
Mo$cibrodzka et al. 2017). These polarimetric measurements 
allow us to carry out new quantitative tests of horizon-scale 
scenarios for accretion and jet launching around the 
М87* black hole. In this Letter we present our interpretation 
of the EHTC VII resolved polarimetric images of the ring 
in M87”. 

Our analysis is presented as follows. In Section 2 we report 
polarimetric constraints from M87* ЕНТ 2017 and supplemen- 
tal observations, and argue that they can be used for scientific 
interpretation, focusing on several key diagnostics of the degree 
of order and spatial pattern of the polarization map. In 
Section 3 we present one-zone estimates of the properties of 
the synchrotron-emitting plasma. In the transrelativistic para- 
meter regime relevant for the M87 core, a full calculation of 
polarized radiative transfer using a realistic description of the 
plasma properties is essential. To that end, in Section 4 we 
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describe a set of numerical simulations of magnetized accretion 
flows that we use to compare with our set of observables. In 
Section 5 we show that only a small subset of the simulation 
parameter space is consistent with the observables. The favored 
simulations feature dynamically important magnetic fields. We 
discuss limitations of our models in Section 6, and discuss how 
future EHT observations can further constrain the magnetic 
field structure and plasma properties near the supermassive 
black hole’s event horizon in Section 7. 


2. Polarimetric Observations and Their Interpretation 
2.1. Conventions in Observations and Models 


Throughout this Letter we use the following definitions and 
conventions for polarimetric quantities, following the Interna- 
tional Astronomical Union (IAU) definitions of the Stokes 
parameters (Z, О, И, V) (Hamaker & Bregman 1996; 
Smirnov 2011). The complex linear polarization field P is 


P = Q + id. (1) 


Then, the electric-vector position angle (EVPA) is defined as 
EVPA = Sa?) (2) 


The EVPA is measured east of north оп the sky. Therefore, 
positive Q is aligned with the north-south direction and 
negative O with the east-west direction. Positive М is at 
a+45 deg angle with respect to the positive O axis (in the 
northeast-southwest direction). Positive Stokes Y indicates 
right-handed circular polarization, meaning in our convention 
that the electric field vector of the incoming electromagnetic 
wave is rotating counter-clockwise as seen by the observer. In 
the synchrotron radiation models that we consider, a positive 
value of emitted Stokes Y is associated with an angle Өв 
between the wavevector k” and magnetic field b" as measured 
in the frame of the emitting plasma in the range 05 c (0, 0.57]. 
Negative Y corresponds to 05 c (0.57, т]. 

The linear and circular polarization fractions at a point in the 
image are defined as 


— IPI 
|т| = I (3) 
IVI 
= : 4 
|| т (4) 


We also define the rotation measure between two wavelengths 
А, and № 
EVPA(A) — EVPA(A2) 


RM = 
A-A 


(5) 


Unresolved observations measure the net (image-integrated) 
polarization fractions 


(Sey + (5и 
y» | 
YM 

Vnet = Уд , 


(6) 


Imlnet = 


(7) 
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Table 1 
ALMA-only Measurements of M87*'s Unresolved Polarization Properties at 
v = 221 GHz (Goddi et al. 2021) 


Day F [maer [уе RM 

dy) (%) (%) (10? rad m?) 
April 5 1.28 + 0.13 2.42 + .03 «0.2 (0.6 + 0.3) 
April 6 1.31 0.13 2.16 + .03 <0.3 (1.5 + 0.3) 
April 10 1.33 + 0.13 2.73 + .03 «0.3 (—0.2 + 0.2) 
April 11 1.34 + 0.13 2.71 + .03 «0.4 (—0.4 + 0.2) 


where the sums are over all pixels 7 in the resolved image. In 
addition to the signed circular polarization fraction vnet, we also 
frequently consider the absolute value |vnetl, as circular 
polarization measurements of the M87" core at 230 GHz do 
not constrain its sign (Goddi et al. 2021). 

In describing the resolved linear polarization in EHT images, 
we define the image-average linear polarization fraction, 
weighted by the total intensity of each image pixel, as 


259i * i 
(Im) = ESI ED (8) 


Note that (||) depends on the imaging resolution (beam size), 
while |m|net is the usual unresolved linear polarization fraction 
and does not depend on resolution. 


2.2. Unresolved Polarization and Rotation Measure 
Measurements toward M87's Core from ALMA 


As part of the EHT 2017 observation campaign, we obtained 
ALMA array measurements of the unresolved, net, near 
230 GHz, polarimetric properties of M87's core and jet on 
2017 April 5, 6, 10, and 11 (hereafter these observations are 
referred to as ALMA-only observations). ALMA-only mea- 
surements are given at v — 221 GHz, a central frequency of 
ALMA Band 6, which has four spectral windows, each 
centered at 213, 215, 227, and 229 GHz. These results, along 
with details on the observations and data calibration, are 
presented in Goddi et al. (2021); we summarize them here in 
Table 1. From the ALMA-only data, the net linear polarization 
fraction (Equation (6)) of the core is |m|net 2 2.796. The data 
also provide an upper limit on the net circular polarization 
fraction (Equation (7)) of the core of |vjna 0.396, with a 
magnitude and sign that vary over the four observed epochs. 
Goddi et al. (2021) also measured an EVPA that varies with 
wavelength across the ALMA band; the slope of EVPA with 
wavelength is consistent with EVPA ox А2, as expected for 
Faraday rotation. The inferred rotation measure (Equation (5), 
for frequencies/wavelengths in ALMA Band 6) is also time 
variable and changes sign between 2017 April 5 and 11, with a 
maximum magnitude [RM] = 1.5 x 10° rad m 2. 

The ALMA-only measurements include extended ~arcse- 
cond scale structures that are entirely resolved out of the EHT 
maps of M87’s core region. As a result, the total 221 GHz flux 
density of M87* measured by ALMA alone is a factor of ~2 
larger than that captured by the resolved EHT images (see 
also EHTC IV). For that reason, we adopt a more conservative 
upper limit of |V|net < 0.896, which would account for the case 
where the large-scale emission is not circularly polarized. 


EHT Collaboration et al. 


2.3. Spatially Resolved Linear Polarization of M87's Core in 
EHT 2017 Data 


The resolved polarimetric images of the M87 core reported 
in EHTC VII display robust features between different image- 
reconstruction algorithms and across four days of observations 
(2017 April 5, 6, 10, and 11). At 20 yas resolution, those 
images consistently show a region of the highest linear 
polarized intensity in the southwest portion of the ring, with 
a fractional linear polarization |m| < 3096 at its maximum. The 
image-average linear polarization fraction takes on values 
5.7% < (|m|) < 10.7% across the different observation days 
and image-reconstruction techniques. The range of the image- 
integrated net polarization fraction is 1.0% < |m|net < 3.7% 
(see EHTC VII, their Table 2). Because polarized emission 
outside the EHT field of view but inside the ALMA-only core 
is unconstrained, we adopt ће ЕНТ [|m|,4 range when 
evaluating models. 

The EHT images on all four observing days reveal a 
characteristic azimuthal pattern of the EVPA angle around the 
emission ring. To quantify this pattern across the image, we 
use a decomposition of the complex linear polarization 
TP = Q + il into azimuthal modes with complex coefficients 
Өл (Palumbo et al. 2020, hereafter PWP). For a polarization 
field in the image plane given in polar coordinates (p, x), the 
бл coefficients аге 


І р (27 ino 
el Peee"wep 09) 


апп Y Pmin 


where Jann is the Stokes Z flux density contained inside the 
annulus set by the limiting radii Phin and pmax We take 
Pmin = О and Pmax to be large enough to include the entire ЕНТ 
image. 

Within the library of polarized images from general 
relativistic magnetohydrodynamic (GRMHD) simulations pre- 
sented in EHTC V, PWP found that the m= 2 coefficient, (5, 
was the most discriminating in identifying the underlying 
magnetized accretion model. The phase of 3, maps well onto 
the qualitative behavior expected of polarization maps with 
idealized magnetic field configurations. In our convention, 
radial EVPA patterns have positive real (5 (L0, = 0 deg), 
azimuthal EVPA patterns have negative real 0 (40 = 
180 deg), and left- (right-) handed spiral patterns have positive 
(negative) pure imaginary 6z (Z05— 90deg and —90 deg, 
respectively). These idealized EVPA pattern configurations and 
their corresponding б, coefficients are summarized іп 
Appendix A and Figure 1 of PWP. The measured range of 
the complex 6 coefficient across the different image- 
reconstruction methods and observing days reported in EHTC 
VIL their Table 2, is 0.04 < || < 0.07 for the amplitude and 
— 163 deg < arg[85] < —129 deg for the phase. 

Appendix A demonstrates that the information in the (5 
coefficient can be obtained in the visibility domain using 
measurements of the linear polarization E — (gradient) and 
B — (curl) modes of the polarization field (e.g., Kamionkowski 
& Kovetz 2016). Trends in @ metric computed across the 
GRMHD image library (Section 4) can be obtained in the 
visibility domain using only E-and B-mode measurements 
taken on EHT 2017 baselines, as long as the visibilities are 
accurately phase calibrated. Because accurate phase calibration 
of EHT data is non-trivial and requires fully modeling the 
polarized source structure, we use image-domain comparisons 
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Table 2 
Parameter Ranges for the Ouantities Used in Scoring Theoretical Models in this 
Letter 

Parameter Min Мах 
ты в CT 
ШЕ 0 0.8% 
(тј) 5.7% 10.7% 
185 0.04 0.07 
ZB —163 deg —129 deg 


Note. The ranges for |т|һеь (|т|), and 5 were taken from EHTC VII Table 2. 
These ranges represent the minimum — 1с error bound and maximum + 1с 
error bound across five different image-reconstruction methods, and incorpo- 
rate both statistical uncertainty in the polarimetric image reconstruction and 
systematic uncertainty in the assumptions made by different reconstruction 
algorithms. The upper limit on |v|;4 was taken as ~2x the maximum value 
found by Goddi et al. (2021). 


to the reconstructions presented in EHTC VII for the 
constraints in this Letter. 

As in total intensity, both the unresolved and resolved 
polarimetric properties of the 230 GHz M87* image changed 
over the week between 2017 April 5 and April 11. Notably, the 
integrated EVPA in the EHT image changes by 2230-40 deg 
(while the ALMA-only EVPA changes by <10 deg). We will 
not interpret this variability further in this work; however, 
Section 6 discusses expectations for time variability from 
viable simulation models. The observational ranges of the key 
parameters that we use in comparing theoretical models to data 
in Section 5—namely |m|net, |У|пеь (|m]), and 82 amplitude and 
phase—are summarized in Table 2. 


2.4. External and Internal Faraday Rotation 


Faraday rotation in a uniform plasma with rotation measure 
(RM) rotates the EVPA away from its intrinsic value EVPA, 
according to Equation (5). The change in EVPA from its 
intrinsic value at 230 GHz (A œ 1.3 mm) is 


RM 
AEVPA | йб —) ce. (10) 
Polarized light rays passing through a uniform medium are 
subject to the same RM. The net source polarization angle is 
then coherently rotated away from its intrinsic value without 
any depolarization. This scenario of “external” Faraday rotation 
has been used to infer the mass accretion rate for sources where 
an RM is measured or constrained (e.g., Bower et al. 2003; 
Marrone et al. 2006, 2007; Kuo et al. 2014), by assuming that 
the observed radiation passes through the bulk of the accretion 
flow. Because relativistic electrons suppress the Faraday 
rotation coefficient аз œ 1/ T? (e.g., Jones & Hardee 1979), 
these models assume that the RM is produced outside the 
emission region at the radius where Ө, = kT,/ m,c? — ], usually 
r~ 100 ғ, (where rz = GM/ c? is the gravitational radius). 
However, in accreting systems like M87", it is unclear 
whether this external Faraday rotation model is a good 
approximation. As we estimate below, one-zone emission 
models of M87* predict substantial RM within the emission 
region itself at radii r S 5 r,. At its low viewing inclination, the 
observed polarized radiation emitted near the horizon may not 
pass through the bulk of the high-density, infalling gas. 
Therefore, “internal” Faraday rotation, which can depolarize 
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the emission as well as rotate the EVPA (Burn 1966), is also an 
important effect to consider. 

The observed ~10% linear polarization of the ring at the EHT 
scale of ~20 pas is much lower than the intrinsic synchrotron 
polarization fraction 270% expected locally. This could result from 
synchrotron self-absorption of the emitted radiation, but one-zone 
estimates and theoretical models (e.g., EHTC V, and references 
therein) suggest that the 230 GHz emission is mostly optically thin. 
It is more likely that the observed depolarization of the resolved 
emission could be the result of polarization structure that is 
scrambled at resolutions finer than the EHT beam. Turbulent 
magnetic fields and Faraday rotation internal to the emission region 
could produce this scrambling. In Section 4.3 we show that 
turbulence in GRMHD models alone is insufficient to produce this 
level of depolarization. Significant internal Faraday rotation of 
polarization vectors on different rays by different amounts can 
produce a sufficiently scrambled image that is depolarized when 
spatially averaged over a telescope resolution element (beam,e.g., 
Burn 1966; Agol 2000; Quataert & Gruzinov 2000; Beckert & 
Falcke 2002; Ruszkowski & Begelman 2002; Ballantyne et al. 
2007). 

From the simultaneous ALMA-only M87" observations, the 
RM implied by changes in the EVPA across the ALMA band is 
[RM] £ 1.5 x 10 rad m ?. These values are consistent with, 
but much more constraining than, the range determined from 
past SMA observations (—3.4—7 х 10° rad m ?, Kuo et al. 
2014). The ALMA-only EVPA difference varies by order unity 
in magnitude and sign over the observing campaign, and 
includes a large flux contribution from extended emission not 
captured by EHT 2017 imaging (EHTC IV). Using a two- 
component model, Goddi et al. (2021) show that the RM 
toward the core emission in the EHT field of view could exceed 
the RM computed from the ALMA-only data, with allowed 
values as large as [RM|<10°radm ?. Because of this 
uncertainty, we do not use the observed RM as an observa- 
tional constraint in our analysis. We account for uncertainty 
related to the observed time variability by using reconstructed 
polarized EHT images from both 2017 April 5 and 11 to define 
the acceptable ranges (see Table 2) of the observational 
parameters used to score theoretical models in Section 5. 


3. Estimates and Phenomenological Models 


In this Section, we take a first look at the importance of 
internal Faraday rotation and magnetic field structure in 
determining the characteristics of the 230 GHz EHT image. In 
Section 3.1 we obtain order-of-magnitude estimates of the 
plasma properties in M87*by interpreting the observed 
depolarization as entirely due to the effect of internal Faraday 
rotation on small scales. In Section 3.2 we explore the effects of 
different idealized magnetic field configurations on the 
observed polarization pattern from plasma orbiting a black 
hole in the absence of Faraday effects. 


3.1. Parameter Estimates from One-zone Models 


Based on a one-zone isothermal sphere model, EHTC V 
derived order-of-magnitude estimates of the plasma number 
density n, and magnetic field strength В in the emitting region 
around M87" as constrained by the Stokes 7 image’s bright- 
ness, size, and total flux density: 


n, ~ 2.9 x 104 cm, (11) 
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Figure 2. Allowed parameter space in number density and dimensionless electron temperature (ne, Ө,) (red region) for the one-zone model described in Section 3.1 for 
three constant values of 3, = 8тп„т„с”Ө,„ / B?. We require that the optical depth 7; « 1 (green region), the Faraday optical depth тру > 2m (blue region), and the total 
flux density 0.2 « F, « 1.2 Jy (black region). Contours of constant magnetic field strength are denoted by labeled dashed lines. 


В = 49 С. (12) 


In this model, the emission radius was assumed to be r œ 5r,, 
and the electron temperature was assumed to be Т. = 
6.25 x 10? К, based on the observed brightness temperature 
of the EHT image. This temperature corresponds to O,— 
КТ, тс? = 10.5, so the emitting electrons һауе moderately 
relativistic mean Lorentz factors 7 ~ 30, ~ 30. The angle 
between the magnetic field and line of sight is set at 0 = 7/3. 
This model ignores several physical effects that are included in 
more sophisticated models and simulations and which are 
necessary to fully describe the emission from M87*. The 
plasma is considered to be at rest and so there is no Doppler 
(de)boosting of the emitted intensity from relativistic flow 
velocities. Redshift from the gravitational potential of the black 
hole is also not included. 

Given ne, B, and Te, we can estimate the strength of the 
Faraday rotation effect at 230 GHz quantified by the optical 
depth to Faraday rotation т, : 


Ty, Sr X py ~ 5.2 E (13) 
8 


where py is the Faraday rotation coefficient (e.g., Jones & 
Hardee 1979). For emission entirely behind an external 
Faraday screen, 7,, is related to the rotation measure RM via 
То, = 2RMA, which follows from the radiative transfer 
equations for spherical Stokes parameters in the absence of 
other effects (see e.g., Appendix A of Moécibrodzka et al. 
2017) and the fact that py ox А. 

Our estimated 7,, indicates that Faraday rotation internal to 
the emission region is an important effect and could thus 
explain the depolarization observed in M87". Faraday effects 
are even more important for the case of polarized light emitted 
by relativistic electrons that travel through a dense, colder 
accretion flow (e.g., Mo$cibrodzka et al. 2017; Ricarte et al. 
2020). In addition, for the same parameters, Faraday conver- 
sion of linear to circular polarization may also be important 
(т ~ 0.5), while self-absorption is weak (7; 0.05). Requir- 
ing an internal Faraday optical depth т, > 27 (large enough to 
produce significant depolarization) provides an additional 


constraint on one-zone models independent of those used 
in EHTC V, which fixed the electron temperature at an 
assumed value. Assuming 7,, > 27 allows us to break the 
degeneracy between magnetic field strength, electron temper- 
ature, and plasma number density. 

Hence, we consider the same model as in EHTC V at several 
different values of = 8тп„КТ„/В”, constrained Ьу Ше 
requirement that the Faraday optical depth 7, > 27. To be 
consistent with the 230 GHz EHT data, we also require that the 
observed image have a total flux F, between 0.2 and 1.2 Jy, and 
that the model has a maximum total intensity optical depth 
7, « 1. Figure 2 shows what constraints these requirements put 
on the electron number density n, and the dimensionless 
electron temperature Ө, at three different values of 8e. For 
values of 0.01 « 8, « 100, in this simple model the electron 
temperature is constrained to lie in a mildly relativistic regime 
250,5 20 (10? < T, < 1.2 х 10!! К), and the magnetic field 
strength is 1 X B < 30 G. The number density of the emitting 
electrons depends more sensitively on the assumed value of у, 
taking on values between 10* cm ? and 10' cm”. 

The one-zone model estimates suggest that both the total 
intensity and polarized emission can be produced in a mildly 
relativistic plasma in a magnetic field of relatively low strength 
B < 30 С. For higher values of B, the electron temperature 
would be too small to explain the observed maximum 
brightness temperature (~10'° К) in the M87* EHT image 
(EHTC IV). Very high values of B are independently 
disfavored by the small degree of circular polarization 
|у. 1% seen in M87*. For В ~ 100 G, the ratio of the 
Stokes Y emissivity to the Stokes Z emissivity jv], i œ 1%. For 
В > 105 G, jv/ji = 10%, for all temperatures >10 ~ K. We also 
note that magnetic fields of B = 5 G are sufficient to produce jet 
powers of Pj. Z 102 erg E (e.g., EHTC V) via the Blandford 
& Znajek (1977) process. 


3.2. EVPA Pattern and Field Geometry 


To demonstrate how the intrinsic magnetic field structure in 
the emission region influences the observed polarization 
pattern, in this section we present the polarization configura- 
tions from three idealized magnetic field geometries around a 
black hole—a purely toroidal field, a purely radial field, and a 
purely vertical field— as seen by a distant observer. In Figure 3 
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Figure 3. (a) Numerical calculations of the polarization configuration generated by an orbiting emission region in the shape of a torus at 8r, in three imposed magnetic 
field geometries and viewed at i = 163 deg (with material orbiting clockwise on the sky). The orbital angular momentum vector is pointing away from the observer 
and to the east (to the left). Total intensity is shown in the background with higher brightness temperature regions shown as being lighter in color. In the foreground, 
the observed EVPA direction is shown with white ticks, with the tick length proportional to the polarized flux. (b) Analytic calculations of the polarization 
configuration from a thin ring of magnetized fluid at 8r, inclined by 163 deg to the observer in the same magnetic field geometries as in (a). While the distribution of 
emitting material is different in the two models, both the sense of asymmetry in the brightness distributions and the polarization patterns match those from the 
numerical calculations. (c) Schematic cartoons showing the emitting frame wavevector k, magnetic field direction Ê, and polarization vector P = Ё x B for each case. 
In the bottom-right panel, k' denotes the approximate light bending contribution to the wavevector. 


we show polarimetric images from these simple field configura- 
tions computed with two methods: a numerical model of an 
optically thin emission region around the black hole (top row of 
Figure 3), and an analytic treatment of the parallel transport of the 
polarization vector that is originally perpendicular to the magnetic 
field (R. Narayan et al. 2021, in preparation, middle row of 


Figure 3). We show the polarization maps from both methods for 
the three idealized magnetic field configurations viewed at an 
inclination angle of i= 163 deg. Both the analytical and 
numerical calculations assume a zero-spin black hole (Schwarzs- 
child metric), though we have found that the main features of 
these polarization patterns are insensitive to spin. 
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In the top row of Figure 3 we show the result of numerical 
calculations performed with the general relativistic ray tracing 
code grtrans (Dexter & Agol 2009; Dexter 2016) of 
polarized emission from an optically and Faraday-thin compact 
emission region, or “hotspot”, in Keplerian orbit around a black 
hole in the equatorial plane. The hotspot has a radial extent of 
3r, and moves in an imposed and idealized magnetic field 
geometry in a circular orbit at a radius of 8r, (following 
Gravity Collaboration et al. 2018, 2020). We construct a 
phenomenological model of a torus of emitting, rotating plasma 
by studying the time-averaged polarized emission images from 
one revolution of this hotspot around the black hole. We have 
verified that a semi-analytic implementation (Broderick & 
Loeb 2006) of a hot accretion flow model (Yuan et al. 2003) 
produces consistent polarization patterns when using the same 
field geometry. 

In the second row of Figure 3, we compare these numerical 
results to results from an analytic calculation of the observed 
polarization pattern generated by the emission of polarized 
light on a thin ring of radius 8r, in the equatorial plane. In this 
model (R. Narayan et al. 2021, in preparation) the polarization 
vectors are emitted perpendicular to the imposed magnetic field 
geometry in the fluid rest frame; they are transformed on their 
way to the observer using an approximate, analytic treatment of 
the effects of light bending, parallel transport, and Doppler 
beaming. This calculation includes radial inflow as well as 
rotation in the velocity field; the models shown use purely 
toroidal motion (clockwise on the sky) with the same idealized 
magnetic field geometries as in the numerical case. The models 
match the asymmetric brightness distributions and polarization 
patterns of the numerical calculations. In particular, both 
models produce consistent helical EVPA pattern in the case of 
a vertical magnetic field. 

The linear polarization direction P of synchrotron radiation 
in the emitted frame is perpendicular to the wavevector k and 
the magnetic field vector B. We define the toroidal magnetic 
field as consisting only of magnetic field components in the 
azimuthal direction, while the poloidal magnetic field consists 
of the remainder, including both radial and vertical compo- 
nents. In a purely toroidal field case, the EVPA shows a radial 
pattern (left column in Figure 3). Purely radial magnetic fields 
(middle column) give a complementary result; the polarization 
has a toroidal configuration, similar to a 90 deg rotation of the 
linear polarization ticks from the toroidal case. 

In a vertical magnetic field (right column in Figure 3), we 
might expect that P should be vertical (north—south) every- 
where since a vertical В is tilted east—west for this viewing 
geometry. We might also expect that P ~ 0 when the black 
hole is viewed face on, because 418. Instead, the linearly 
polarized emission from a purely vertical field shows a twisting 
pattern that wraps around the black hole. The twist results from 
a combination of light bending and relativistic aberration. Light 
bending in the emitting region near the black hole contributes a 
radial component k' to the emitted wavevector Å that initially 
points away from the black hole (see the schematic cartoon in 
the bottom-right panel of Figure 3). As a result, close to the 
black hole, the total wavevector kemi = k +k’ and the 
magnetic field B are no longer parallel, the polarization is 
non-zero, and the resulting EVPA pattern is north—south 
symmetric. Relativistic motion of the emitting material 
(aberration) breaks the symmetry and gives the twisting pattern 
a handedness corresponding to the orbital direction. For the 


EHT Collaboration et al. 


pure vertical field considered here, the handedness depends on 
the rotation direction and the observed pattern is consistent 
with clockwise rotation. The dependence on direction of 
motion and magnetic field configuration are discussed in more 
detail in a forthcoming paper (R. Narayan et al. 2021, in 
preparation). The EVPA patterns in these images do not show a 
strong dependence on the black hole spin. 

In a rotating flow, weak magnetic fields are sheared into a 
predominantly toroidal configuration (e.g., Hirose et al. 2004). 
In the absence of other effects (e.g., external Faraday rotation), 
the observed azimuthal EVPA pattern suggests the presence of 
dynamically important magnetic fields in the emission region, 
which can retain a significant poloidal component in the 
presence of rotation. In the next sections, we compare 
numerical simulations of the accretion flow and jet-launching 
region in M87" with different field configurations to the 
EHT2017 data to better constrain the magnetic field structure. 


4. M87* Model Images from GRMHD Simulations 


The low resolved fractional linear polarization observed by 
the EHT contradicts the results from an idealized magnetic field 
structure with no disorder. For typical parameters of the 
230 GHz emission region, Faraday rotation and conversion are 
expected to be important. Magnetic field structure, plasma 
dynamics and turbulence, and radiative transfer effects 
including Faraday rotation can be realized in images from 
three-dimensional general relativistic magnetohydrodynamic 
(3D GRMHD) simulations of magnetized accretion flows. We 
use 3D GRMHD simulations (described in Section 4.1) in 
combination with polarized general relativistic radiative 
transfer (GRRT) models (described in Section 4.2) to model 
polarized images of M87". In Section 4.3, we describe trends of 
the key observables (||, |у|һеь (|m|), and (2) in our 
GRMHD polarimetric image library. 


4.1. GRMHD Model Description 


The simulation library generated for the analysis of the 
EHT 2017 total intensity data in EHTC V consists of a set of 
3D GRMHD simulations that were postprocessed to generate 
simulated black hole images via GRRT. For simulations using 
black holes with non-zero angular momentum, we only 
considered accretion flows in which the angular momentum 
of the flow and the hole were aligned (parallel or anti-parallel). 
Because the equations of non-radiating ^? GRMHD are scale 
invariant, each fluid simulation was thus fully parameterized by 
two values describing the angular momentum of the black hole 
and the relative importance of the magnetic flux near the 
horizon of the accretion system. A comparison of several 
contemporary GRMHD codes, including those used to generate 
the simulation library, can be found in Porth et al. (2019) and in 
H. Olivares et al. (2021, in preparation). 

The black hole angular momentum J is expressed in terms of 
the dimensionless black hole spin parameter a, = Jc/ GM. In 
this Letter, we consider simulations run with the iharm code 
(Gammie et al. 2003; Noble et al. 2006) with a, — —0.94, 
—0.5, 0, 0.5, and 0.94, where positive (negative) spin implies 
alignment (anti-alignment) between the accretion disk and the 
black hole angular momentum. Several studies of “tilted” disks 


129 We assume that M87* can be described by models in which radiative 


cooling is negligible so that it does not affect the dynamics of the plasma and 
images can be generated in postprocessing. 
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have been conducted (e.g., Fragile et al. 2007; McKinney et al. 
2013; Morales Teixeira et al. 2014; Liska et al. 2018; White 
et al. 2019; Chatterjee et al. 2020). As there does not yet exist a 
full library of tilted disk simulations spanning a range of spins, 
we limit our analysis to the aligned and anti-aligned 
simulations considered in EHTC V. 

The strength of the magnetic flux near the horizon 
qualitatively divides accretion flow solutions into two cate- 
gories: the Magnetically Arrested Disk (MAD) state (e.g., 
Bisnovatyi-Kogan & Ruzmaikin 1974; Igumenshchev et al. 
2003; Narayan et al. 2003) in which the magnetic flux near the 
horizon saturates and significantly affects the dynamics of the 
flow, and the contrasting Standard and Normal Evolution 
(SANE) state (e.g., De Villiers et al. 2003; Gammie et al. 2003; 
Narayan et al. 2012). The relative importance of magnetic flux 
in a simulation is quantitatively described by the dimensionless 
quantity 


ф = égg(Mr2c) 1/2, (14) 


where Фвн is the magnitude of the magnetic flux crossing one 
hemisphere of the event horizon (see Tchekhovskoy et al. 
2011; Porth et al. 2019) and M is the mass accretion rate 
through the event horizon. The flux saturates at values of 
ф 2, 50, and the flow becomes MAD. The SANE simulations 
that we consider have lower values of фә 5,90 Accreted 
material supplied at large scales could, in principle, supply any 
value of net vertical flux. Here, we do not explore cases with 
small or zero net vertical flux ф < 1. We also do not consider 
values in the relatively narrow intermediate range 5 < ф < 50. 

The SANE simulations considered here used a grid 
resolution of 288 x 128 x 128, a fluid adiabatic index y= 4/ 
3, and an outer simulation domain of Fow = 50 r. The MAD 
simulations used a grid resolution of 384 x 192 x 192, an 
adiabatic index ^ = 13/9, and an outer simulation domain of 
Tout 10? rg. The simulations were carried out in modified 
spherical polar Kerr-Schild coordinates, where grid resolution 
is concentrated toward the midplane to help resolve the 
magnetorotational instability. АП models in the ЕНТ library are 
evolved for at least г = 10* r,/c and their accretion flows reach 
steady state within r= 10-20 r,. 


4.2. Ray-traced Polarimetric Images from GRMHD 
Simulations 


Unlike the equations of GRMHD, the equations of radiative 
transfer are not scale invariant, and so we must introduce two 
scales to the simulation when we ray-trace images from the 
numerical fluid data. The length (and time) scale is set „БУ Ше 
mass of the black hole, assumed to be Mgy = 6.2 х 10? Mo in 
accordance with the value used to generate the EHTC V 
simulation library. For our models, we also adopt the D — 16.9 
Mpc distance to M87* used in EHTC V. The density scale of 
the accreting plasma (equal to the scale of the magnetic 
pressure) is chosen so that on average the simulated images 
pas the observed 230GHz compact flux density, 

F,c 0.5 Ју. 

Images were generated from the set of simulations for 

several values of the polar inclination angle i chosen to be 


130 Note that the MAD threshold Oz 2 50 is given in Gaussian units where 


[Ф] = Gem”. If the field strength is given in the Lorentz—Heaviside units 
typically used in simulations (Big = BG / X47), the MAD threshold on the 
dimensionless flux is ф œ 15. 
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broadly consistent with observational estimates of the inclina- 
tion angle of the M87 jet (e.g., Walker et al. 2018). The 
position angle on the sky can be changed after image 
generation by rotating both the image and the Stokes O and 
U components appropriately. Each image has a 320 x 320 
pixel resolution over a 160 pas field of view, where each pixel 
contains full Stokes Z, Q, U, V intensities. 

In GRMHD simulations, we make the approximation that the 
plasma is thermal, i.e., that the electrons and ions are described 
by a Maxwell-Jüttner distribution function (Jüttner 1911). 
However, the plasma around M87* and in other hot accretion 
flows is most likely collisionless, with electrons and protons 
that are unable to equilibrate their temperatures (e.g., Shapiro 
et al. 1976; Ichimaru 1977). We mimic collisionless plasma 
properties in producing images from the GRMHD simulations 
by allowing the electron temperature 7, to deviate from the 
proton temperature 7;. The simulations used in this work only 
track the total internal energy density gas, not the distinct 
electron and ion temperatures. We set Т. after running the 
simulation according to local plasma parameters following the 
parameterization introduced by Mo$cibrodzka et al. (2016; see 
also Moscibrodzka et al. 2017 and EHTC V). The ratio 
between the ion and electron temperatures R is determined by 
the local plasma p= Pgas/Pmag» where ра = Cy = l)uga,, and 
Pmag = B z 87. The temperature ratio is then taken to be 

2 
R= 2 = Rhigh БЕ 08 + Коу 


- ET (15) 


1-8 
where Ryisn (Riow) are the free parameters of the model and give 
the approximately constant temperature ratio at high (low) 5. 
This approach allows us to associate the electron heating with 
magnetic properties of the plasma. 

In calculating the electron temperature, we further assume 
that the plasma is purely ionized hydrogen and that ions are 
nonrelativistic with an adiabatic index y, = 5/3 and electrons 
are relativistic with у= 4/3. Then, given Ugas from the 
simulation and R from Equation (15), (EHTC V): 

2 MpUgas 


Ti, ЫР ЕВ. (16) 
3pk(2 + В) 


We note that this procedure is not entirely self-consistent, as ће у 
of the combined electron-ion fluid will change depending on the 
relative pressure contributions of electrons and protons while we 
assume it is fixed throughout the simulation domain. See 
Sadowski et al. (2017) for an alternative, self-consistent approach. 

In this Letter, we consider a library of 72,000 simulated 
images composed of sets of 200 realizations of the same 
accretion system described by a fixed set of heating and 
observation parameters. Each set of 200 images is drawn from 
output files spaced by 25—50 r,/c from the set of 10 GRMHD 
simulations spanning five spin values in both MAD and SANE 
field configurations. The inclination angle for each image is set 
to one of either i = 12, 17, 22 deg (retrograde models, а < 0) 
or i= 158, 163, 168 deg (prograde models, a, > 0), according 
to which parity is required to orient the brightest portion of the 
ring in the southern part of the image while ensuring the 
orientation of the approaching jet is consistent with large-scale 
observations. 

We use electron heating parameters Aj, = 1, 10 and 
Ryign — 1, 10, 20, 40, 80, or 160 in Equation (15). EHTC V 
only considered models with Rj,,, = 1. Larger values of А, 
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correspond to lower electron-to-proton temperature ratios in the 
low 2 regions (e.g., the jet funnel). This choice is physically 
motivated for M87", where radiative cooling of the electrons 
may keep T, « T; even in magnetized regions where electron 
heating is efficient (e.g., Mo$cibrodzka et al. 2011; Ryan et al. 
2018; Chael et al. 2019). Lower electron temperatures in 
Коу = 10 models increase the Faraday rotation depth and can 
result in increased depolarization in parts of the image. 

GRMHD simulations produce a highly magnetized jet funnel 
above the black hole's poles, away from the accretion disk. In 
the funnel, where the plasma magnetization parameter 
с=В” / 4трс” > 1, our numerical methods typically fail to 
accurately evolve the plasma internal energy. In the image 
library, we cut off all emission in regions where c > | to ensure 
that we limit the emitting region to plasma whose internal 
energy is safely evolved without numerical artifacts (as 
in EHTC V). We tested the importance of a с> 1 electron 
population by generating a supplementary set of images from 
all models with a cut at c — 10 and found that it did not change 
the overall distribution of the derived metrics we use for model 
scoring in Section 5. 

Each set of 200 model images with the same parameters in 
the image library requires a unique density scaling factor that is 
determined by matching the average flux density from the 
model to the observed compact flux density of M87" measured 
by the EHT. Hence, the mass accretion rates, radiative 
efficiencies, and jet powers will differ between two models 
even if they are derived from the same underlying simulation 
(e.g. if Аһ, Ком, Or i are changed). The additional models 
discussed in Section 6, which explore the effects of different c 
cutoff values and the inclusion nonthermal electrons, also 
require unique mass-scaling factors. 

All of the polarimetric images from GRMHD simulations 
that we analyze in this Letter were generated using the ipole 
code (Moscibrodzka & Gammie 2018), which has been tested 
against analytic solutions and numerical ones produced by 
other numerical GRRT codes (Dexter 2016; Moécibrodzka 
2020). A more comprehensive comparison of various GRRT 
codes that perform total intensity transport and fully polarized 
GRRT can be found in Gold et al. (2020) and B. Prather et al. 
(2021, in preparation), respectively. Preliminary results from B. 
Prather et al. (2021, in preparation) show that the tested codes 
are consistent at the fraction of 1% in all Stokes parameters. All 
calculated images in this work ignore light travel time delays 
through the emission region (the so-called “fast light” 
approach), and are calculated at a single frequency of 
v = 230 GHz, neglecting the finite observing bandwidth of 
the EHT. We confirm that neither of those effects are important 
for models of interest for M87”. 


4.3. Sample GRMHD Model Images and Polarization Maps 


In Figures 4 and 5 we show images and polarization maps 
for a subset of library models. In general, because the horizon- 
scale magnetic fields in MAD models are strong enough not to 
be advected with the accretion flow, they are more likely to 
have a significant poloidal component and produce azimuthal 
EVPA patterns (Figure 3). In contrast, SANE models tend to 
show more radial EVPA patterns. Some MAD a, = 0.94 and 
SANE a, = 0 images are notable exceptions to this trend. 
These trends are also apparent in the distributions of the (5 
phase across the full image library that we consider later in 
Figure 9. 
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The GRMHD models at their native resolution include 
notable disorder in the EVPA structure, resulting from both 
magnetic turbulence and Faraday rotation. Models with larger 
Ryign have lower electron temperatures and higher Faraday 
rotation depths, resulting in the most disordered polarization 
maps. Many of the EVPA patterns seen in the images blurred 
with a 20 uas Gaussian kernel to simulate the limited ЕНТ 
resolution resemble those from the idealized magnetic field 
models in Figure 3, indicating that the net EVPA pattern after 
blurring may trace the intrinsic magnetic field structure. 

In Figure 6 we show a sample polarization map at full 
resolution compared to the same map blurred with circular 
Gaussian kernels of 10 pas and 20 uas FWHM. From tests with 
synthetic data, blurring (convolving with a circular Gaussian 
kernel) provides a reasonable approximation to image recon- 
struction from the EHT data at a comparable resolution (EHTC 
VID. The resolved average fractional polarization in the blurred 
images (|m|) traces the degree of order in the intrinsic 
polarization map. In the blurred images, disordered polarized 
structure on small scales produces beam depolarization. The 
degree of depolarization decreases with increasing spatial 
resolution (decreasing beam size). 

The bottom row of Figure 6 shows the same unblurred and 
blurred polarization maps, but calculated without the effect of 
Faraday rotation (py — 0). Those images show more coherent 
EVPA structure, with much larger |m|;« and, particularly when 
blurred, much larger (|m|). Evidently, for this particular model, 
the depolarization visible in the corresponding top panels is due 
to Faraday rotation internal to the emission region. In addition, 
the net EVPA pattern shifts by a significant amount. The 
change in 62 by ~80 deg would correspond to an apparent RM 
of ~—4 x I0? rad m ?. Our GRRT calculations include all 
Faraday rotation occurring inside the GRMHD simulation 
domain (rout = 50-100 r,), both external and internal to the 
230 GHz emission region. The observables considered here, for 
the low viewing inclination of M87*, do not depend strongly on 
that outer radius, as long as it is at r 2 40 r,. We cannot rule out 
the presence of additional Faraday rotating material at larger 
radii 2100 r,, and its effects are not included in our models. 
Appendix B discusses the origin of the RM in our models in 
more detail. 


4.4. GRMHD Model Theory Metrics 


We compute the polarimetric observables (||, |V|net: 
(|m|), G2) described in Section 2.3 from model images blurred 
with a circular Gaussian kernel with a FWHM of 20 was in 
order to compare them to the ranges measured from EHT and 
ALMA-only data. Both (|m|) and (2 depend on the resolution 
and hence the size of the Gaussian blurring kernel. The value of 
б» also depends on the choice of the image center. We do not 
shift the library images before computing Bm coefficients for 
comparison with the range inferred from the EHT image 
reconstructions, which have been centered by aligning them to 
the centered, fiducial total intensity images in EHTC IV. As 
discussed in Palumbo et al. (2020), a centering offset u 
expressed as a fraction of the diameter of a PWP m= 2 ring 
causes a quadratic falloff in 32 power as 662/|62| ez 4u’. Effects 
on the (5 phase enter at similar order. In the case of ће ЕНТ 
image, u is likely less than one-fifth, meaning that centering 
errors in (25 will be sub-dominant to other uncertainties, such as 
the choice of the blurring kernel or the variation across methods 
and days. 
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Figure 4. Sample snapshot false-color images and polarization maps for a subset of the models in the EHT M87* simulation image library at their native resolution 
(top three rows) and blurred with a 20 pas circular Gaussian beam (bottom three rows). The inclination angle for all images is either 17 deg (for negative a, models) or 
163 deg (for positive а, model), with the black hole spin vector pointing to the left and away from the observer. The tick length is proportional to the polarized flux, 
saturated at 0.5 of the maximum value in each panel. Here models with Ку, = 1 are shown. In general, the EVPA pattern is predominantly azimuthal for MAD 


models (e.g., MAD a, = 0 Rnign = 1) and radial for SANE models (e.g., SANE а, = 0.94 


= 1), although the SANE a, = 0 models in particular are exceptions 


to this trend. All models show scrambling in the polarization structure on small scales from internal Faraday rotation, with more pronounced scrambling in models 


with cooler electrons (larger Кмьһ parameter). 


Figure 7 (right panel) shows the resolved average polariza- 
tion fraction (|т|) as a function of their image-averaged 
Faraday rotation depth, (7,,). At small (7,,), the average 
polarization fraction is (|m|) = 2096—5096. Intrinsic disorder in 
the magnetic field structure due to turbulence is generally 
insufficient to produce the low observed image-average 
polarization fraction іп EHT 2017 M87" data (5.796 «(|m|) < 
10.796). This is especially evident for the SANE models with 


prograde black hole spin, which have the highest resolved 
polarization fractions. At large (7,,), strong scrambling from 
internal Faraday rotation typically results in small predicted 
polarization fractions of <5% at the scale of the ЕНТ beam. 
The clear exceptions to this trend are some SANE retrograde 
models (a, = —0.9375 for large Rnign), which show (|m|) = 
10%-20% despite their large (т,,) 2, 103. In these models, 
most of the observed polarized flux originates in the forward 
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Figure 5. Same as in Figure 4 but for Rj,,, = 10. We find similar trends, but with more scrambling from larger Faraday depths due to lower electron temperatures. 


jet, while most of the computed Faraday depth is accumulated 
near the midplane. Photons that travel from the forward jet to 
the observer do not encounter the large Faraday depth. For 
similar reasons, the inferred RM can be much lower than 
implied by their large values of integrated т,,. 

Distributions of all observables are shown in Figure 7 ((|mJ), 
left panel), Figure 8 (\m|net and |v|net), and Figure 9 (|| and 
Z5) SANE models tend to have a lower integrated 
polarization fraction and larger circular polarization fraction 
than М87* at 230 GHz. In many cases this is a result of very 
large Faraday rotation internal to the emission region. MAD 
models tend to have larger net linear polarization fraction than 
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observed in М87*. The resolved average fractional polarization 
produces similar trends. Most SANE models with prograde 
spin are too scrambled and most MAD models are too ordered 
compared to the reconstructed polarization maps of M87”. Full 
distributions for all models, including their Ryisn, Riow, and ax 
dependence, are further discussed in Appendix C. 


5. Model Evaluation 
5.1. Model Constraints from Polarimetry 


To evaluate whether a given GRMHD model is consistent 
with the EHT observations reported in EHTC VII, we require 
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Figure 6. Top-left panel: a sample image library polarization map at original resolution, taken from the MAD a, = 0.5 (Riow = 10, Rnigh = 80) model. Top-middle 
and top-right panels: the same map but convolved with a 10 pas and 20 pas FWHM circular Gaussian beam, respectively. The position angle (PA) of the black hole 
spin in all frames is PA = 4-90 deg and the inclination angle isi — 158 deg, meaning that the black hole spin points left and away from the observer. The bottom row 
shows the same model but calculated with py = 0 (no Faraday rotation). When Faraday rotation is excluded, the EVPA pattern is more coherent, resulting in much 
larger values of |т|һе and (|m|). There is also a net rotation of ће EVPA pattern between the two cases, by ~80 deg in the phase of (5. 
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Figure 7. Left panel: distribution of image-averaged fractional polarization (|m|) over the M87* library images blurred with a 20 pas beam. The measured range from 
reconstructed polarimetric images of M87* is shown in dashed lines. Right panel: (|m|) as a function of the intensity-weighted Faraday depth across each image for 
library images blurred with the same 20 pas circular Gaussian beam. The Faraday depth is calculated as the intensity-weighted sum of |py| integrated along each ray. 
For clarity, we show the median (points) and standard deviations (error bars) of the full distributions. The Faraday depth increases monotonically with increasing Rpign 
for fixed values of the other parameters. A large Faraday depth corresponds to scrambling of the polarization map, which decreases the coherence length of the EVPA 
(Jiménez-Rosales & Dexter 2018). Increased scrambling results in stronger depolarization at the scale of ће ЕНТ beam and lower values of (|m|). 


images from the model to satisfy constraints on the four parameters 
derived from the reconstructed EHT images and ALMA-only 
measurements presented in Table 2 and summarized again here. 


2. The image-integrated net circular polarization |v|ne satisfies 
an upper limit from ALMA-only measurements reported in 
Goddi et al. (2021): |v|net < 0.8%. 


1. The image-integrated net linear polarization |m|;« is in 3. The image-averaged linear polarization (|m|) is in the 
the measured range from the EHT image reconstructions: measured range from the EHT image reconstructions at 
1% < |m|y« < 3.7%. 20 pas scale resolution: 5.7% < |m|net < 10.7%. 
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Figure 8. Distributions of image-integrated net linear (left panel) and circular (right panel) polarization fractions for all ЕНТ M87* library images. The dashed lines 
show the allowed range inferred from ЕНТ image reconstructions (for |т|е) and ALMA-only data (for |v|net). 
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Figure 9. Distributions of 3, amplitude (left panel) and phase (right panel) for ЕНТ M87* library images blurred with a 20 pas beam. The measured ranges from 


reconstructed images of M87* are shown as dashed lines. 


4. The amplitude and phase of the complex (5 coefficient 
quantifying coherent azimuthal structure are in the 
measured range: 0.04 < |62| < 0.07 and —163 deg < 
arg[2] < —129 deg. 


We use 72,000 library images (from Section 4) with 200 
time snapshots per model at three inclination angles, six values 
of Аһ = 1, 10, 20, 40, 80, 160, two values of Row = 1, 10, 
five values of а, = —0.9375, —0.5, 0, +0.5, +0.9375, and 
realized with both MAD and SANE magnetic field 
configurations. 

In comparing models to observables, the (2 metric is the 
most constraining. Only 790 snapshot images out of the 72,000 
considered fall in the range of those reconstructed in both (5 
amplitude and phase, compared to 11,526 snapshots for both 
|т|һе and |v|net and 7,727 for the resolved image-average linear 
polarization fraction (|m)). 

Below we explore two quantitative methods for scoring 
models, either by requiring that at least one single snapshot 
image from a model simultaneously passes all constraints 
("simultaneous scoring;” see Section 5.2) or that each 
observational constraint is satisfied by at least one snapshot 
image from a given model (“joint scoring;" see Section 5.3). 
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5.2. Simultaneous Snapshot Model Scoring 


In the simultaneous scoring procedure, we rule out models 
where none of the 600 snapshot images (200 time samples at 
three inclination angles) can simultaneously satisfy the 
constraints on all of the polarimetric observables. Only 73 
out of 72,000 snapshot images across 15 of 120 models 
simultaneously pass all of the constraints. Of those, all but two 
viable snapshot images come from a MAD model. The only 
models with more than five passing images are MAD a, — 0 
Riow =1 Rhign = 160 and MAD ay 0.5 Riow = 1 Ryign = 
80, 160. 

Figure 10 shows three viable snapshot images from both 
SANE and MAD models as well as three snapshot images from 
models ruled out by simultaneous scoring (ie., with no 
snapshots in the entire sample from the model simultaneously 
satisfying all constraints). These images are representative of 
the snapshots that simultaneously satisfy all constraints on the 
image-integrated metrics; they have not been selected based on 
detailed matching of the resolved polarization structure to the 
EHT images. Nonetheless, the top row of images show good 
qualitative agreement with the primary features of the EHT 
image in Figure 1. In contrast, the snapshots from the ruled-out 
models tend to be too polarized, too depolarized, or too radial 
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Figure 10. Sample 230 GHz image library polarization maps shown in the same style as the ЕНТ image in Figure 1. All maps are blurred with a 20 pas circular 
Gaussian beam. In all images, the inclination angle is either 17 deg (for negative a, models) or 163 deg (for positive a, model), with the black hole spin vector 
pointing to the left and away from the observer. The top row displays SANE (a, = 0, Rnign = 80) and two MAD snapshots (both a, = —0.5 and А = 160) from 
left to right. All of the top row images satisfy simultaneous constraints on image-integrated metrics (||, |У|пеь (171), |82|, 22) derived from the EHT2017 results. 
The bottom row displays snapshots from models that fail to produce any images that simultaneously satisfy the observational constraints. These snapshots are from 
two SANE (a, = 0.5 and А. = 1 and 160) and one MAD (a, = 0.94, Къ; = 160) model, from left to right. The failing images are either more polarized than the 
data (left and right panels) or too depolarized (middle panel). They also fail to match the observed EVPA pattern (82 phase). 


in their EVPA pattern. Figure 11 shows the distributions of |82| with x. n E daa" The joint likelihood of each model is the 
for all 600 snapshots from these three passing and three failing product £ = I;£; of those for the five metrics ху. 
models. Although variability is present, the systematic To produce a non-zero likelihood £ in this method, at least 
differences over the five observables considered allow us to one snapshot from a model must lie further from its mean than 
constrain the models. The left panel of Figure 12 shows the the data value does. That can be a different snapshot for each 
total number of images that pass simultaneous scoring as a metric, which makes this method more lenient than the 
function of model, summing over the six Rhign values. simultaneous scoring method. We also note that snapshots 
А СРЕ | are allowed to have the wrong sign ОЁ the difference with their 
5.3. Joint Distribution Model Scoring mean, due to the definition of x* and our use of the mean of the 
In the joint scoring procedure, we use the measured model snapshots themselves. In practice, this makes little 
distributions of the data metrics to ask whether the observed difference in the results. 
value of each metric for M87* is consistent with being drawn In this method, we consider models viable whose joint 
from the distribution seen in the GRMHD simulations. likelihood is 7196 of the maximum found from any model. The 
To do this, we measure X^ values for the five metrics right panel of Figure 12 shows the resulting joint likelihoods 
xj € {|| [уе (Im), [82]. Z 82) for all snapshots k from summed over Rhign- 
a given model as 
5 (худ — ¥)? Ра 5.4. Сотрагіѕоп оў Scoring Results 
jk o^ : qn The results of both scoring procedures are summarized in 
Figure 12, summed over Rpign. Both scoring methods prefer 
where x; , are the values of a scoring metric x; for each of the MAD models to SANE models, with most of the passing 
600 snapshots k from a given model, x; is the mean of those кагы coming from the МАР а„=0 and а, =+0.5 
simulations. 


values for the model, and o; is taken as one half of the observed 
data range from Table 2. Note that the scoring results of this 
method do not depend on the choice of су. We then calculate an 


The main difference between the two scoring procedures is 
that joint scoring prefers Riow = 10 models, while Riow = 1 is 
preferred by simultaneous scoring. SANE models with 


analogous x data Value for the midpoint of the measured range а, = 0.94, Ry, — 1, 10, and Ry, — 10 are ruled ош by 
from Table 2. A likelihood value £; of the data being drawn simultaneous scoring, but score fairly well in joint scoring. For 
from the model distribution is defined as the fraction of images the favored MAD models, when Riow= 1, there are more 
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Figure 11. Distributions of |82| for the sample passing and failing models in Figure 10. The dashed black lines mark the measured values for the snapshot images in 
Figure 10, and the blue bands show the range inferred from ЕНТ M87* data. The models can be constrained using EHT observables even in the presence of significant 
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Figure 12. Results of the simultaneous (left panel) and joint (right panel) scoring methods for comparing GRMHD models to М87* observables. The simultaneous 
scoring method shows the total number of viable images for each image library model after summing over Rhign. Out of a total of 73 passing images, only two are from 


a SANE model. The right panel shows the joint likelihood of each library model after summing over Rpign. In this method, Riow 


SANE a, = +0.94 Rhign = 10 models are also allowed. 


images that simultaneously satisfy all constraints, but when 
Коу = 10, the distributions generally stay closer to the 
observed data ranges and are thus favored by the joint scoring 
method. Due to differences between simultaneous and joint 
scoring results, as well as other methods we have tried, we 
consider the inferred parameters of Riow, Rnigh, and a from 
passing models to be less robust than the overall trend that 
MAD models are favored. 

The simultaneous scoring method has the advantage of 
conceptual simplicity, and of applying each constraint 
simultaneously per image (thus accounting for correlations 
between the scoring metrics). Simultaneous scoring is more 
strict and rules out more models than joint scoring, but it may 
be more limited by the finite number of images generated per 
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10 MAD models are preferred and 


model. The joint scoring procedure has the advantage of being 
more conservative in disfavoring models, but assumes the 
observational constraints are independent in calculating a joint 
likelihood. Instead, they are correlated (in particular |m|net, 
(Iml), and |63). 

The number of images in each model that pass each 
constraint individually (used in joint scoring) and that 
simultaneously pass all constraints (used in simultaneous 
scoring) are presented in Appendix D. 


5.5. Combined EHTC V and Current Polarimetric Constraints 


EHTC V presented constraints on the GRMHD simulation 
models based on fits to the EHT total intensity data, model 
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self-consistency (requiring а radiative efficiency less than that 
of a thin accretion disk at the same black hole spin), and M87's 
measured jet power (requiring a simulation to produce a jet 
power consistent with a conservative lower limit of that from 
M87*, S10” ergs !). Those constraints ruled out MAD 
d, = —0.94 models (from failing to satisfy the EHT image 
morphology), SANE models with a, — —0.5, and all models 
with a, = 0 (from failing to produce enough jet power). Here 
we retain only the jet power lower limit, which is the most 
constraining and straightforward to apply to the expanded 
image library considered in this work. 

Relativistic jets launched in GRMHD simulations (defined 
here as in EHTC V, with a cutoff of Gy > 1) are fully consistent 
with being produced via the Blandford—Znajek process (e.g., 
McKinney & Gammie 2004; McKinney 2006). As a result, 
ағ = 0 models have small or zero jet power, Pje, and are 
rejected by this constraint. These models can still produce 
significant total outflow powers (Pout in EHTC V) in a mildly 
relativistic jet or wind. Many other models with low values of 
Ryign Or moderate black hole spin, including those of SANES 
with a, — 4- 0.94, which are allowed by the joint scoring 
procedure, are also ruled out by the jet power constraint (see 
Table 3 in Appendix D and EHTC V). Combining the 
simultaneous scoring polarimetric constraints with the jet 
power constraint results in 15 remaining viable models: all 
МАР», and spanning the full range of non-zero a, explored. 
This conclusion does not depend on the choice of the 
simultaneous or joint model-scoring procedure. 


6. Discussion 


The resolved EHT2017 linear polarization map of 
M87* shows a predominantly azimuthal linear polarization 
(EVPA) pattern, and relatively low fractional polarization of 
<20% on 20 uas scales. We interpret the low fractional 
polarization as the result of Faraday rotation internal to the 
emission region, which acts to rotate, scramble, and depolarize 
the resolved polarized emission. Adopting this constraint in a 
one-zone model, we estimate typical values of particle density 
ne, magnetic field strength B, and electron temperature Te. In 
semi-analytic emission models with externally imposed, 
idealized magnetic field configurations, azimuthally dominated 
EVPA patterns are produced by poloidal (radial and/or 
vertical) magnetic field components. To fully capture the 
complicated combined effects from magnetic field structure, 
turbulence, relativity, and Faraday rotation on polarimetric 
images of M87", we turn to radiative transfer calculations from 
GRMHD simulations. 

We compared a large image library of emission models from 
GRMHD simulations with metrics designed to capture these 
salient features of the data. The combined constraints of a 
predominantly azimuthal EVPA pattern and a low but non-zero 
fractional polarization are inconsistent with most SANE 
GRMHD models with weaker horizon-scale magnetic fields. 
Some MAD models with relatively cold electrons, realized in 
our library by larger values of Rhigh and/or Куу, remain 
consistent with the data. Here we discuss the implications of 
our results, and limitations in our set of theoretical models that 
may impact our interpretation. 
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6.1. Near-horizon Plasma and Magnetic Field Properties in 
Passing Models 


Both our one-zone and GRMHD models find similar plasma 
conditions in the 230 GHz emission region, driven by the 
requirements of weak 230 GHz absorption and strong 230 GHz 
Faraday rotation. In viable GRMHD models, we find average, 
intensity-weighted plasma properties in the emission region of 
ne~ 10? cm 2, Bc 7-30 С, and 6,~ 8-60. These are in 
good agreement with our one-zone estimates (Section 3.1). We 
have also calculated the intensity-weighted values of the 
absorption and Faraday optical depth, т; and 7,,, over 
snapshots that simultaneously satisfy all our observational 
constraints. The median values are 7; 20.1 and 75, = 50. All 
of our viable images have 7,, > 27, while 2 out of 73 have 
т; 2, 1, consistent with our assumptions in Section 3.1 that the 
plasma Faraday depth is large while the Stokes 7 optical depth 
is small. 

By quantitatively evaluating a large library of images based 
on GRMHD models (Section 5), we identify 25 out of 120 
models that remain viable after applying constraints based only 
on EHT and ALMA-only polarimetric observations. Addition- 
ally applying a cut on jet power of Pje > 102 erg s! (EHTC V) 
rules out the five viable SANE models and all a, = 0 models. 
The precise number and identity of the viable models depends 
mildly on the chosen scoring procedure and on the Gaussian 
blurring kernel size applied to the EHT image reconstructions 
and library simulated images. The overall preference for MAD 
over SANE models is found from both the simultaneous and 
joint scoring procedures, as well as other variants. After applying 
the jet power constraint, no viable SANE models remain for any 
of the scoring methods that we explored. 

MAD models are associated with dynamically important 
magnetic fields. The significant poloidal components of those 
fields can produce a predominantly azimuthal polarization 
pattern (Figure 4), similar to those seen in idealized models 
with prescribed poloidal magnetic fields (Figure 3). Strong 
Faraday effects complicate a direct interpretation of the 
observed EHT polarization map in terms of those idealized 
models. Still, our more detailed comparison favoring MADs 
suggests the presence of dynamically important magnetic fields 
in the emission region on event-horizon scales. 

In Figure 13 we present mass accretion rate and jet power 
distributions both for the viable models identified in EHTC V 
and when adopting the new constraints from polarimetry. ^! 
Polarimetric constraints break degeneracies present in the 
single epoch total intensity data, allowing us to estimate a mass 
accretion rate onto the black hole of M œ (3-20) x 
10-^M., yr-!. This corresponds to й = M/Mgag œ (2-15) x 
10-5, where Мьша is the Eddington accretion rate.?? The 
measured radiative efficiency є = L/Mc? (where L is the 
bolometric luminosity) of the passing models is relatively high 
for a hot accretion flow model: є < 1%. These models have jet 
powers of Pje = 1022743 erg s™!. 

The mass accretion rate found here is much lower than 
the Bondi rate calculated from Chandra observations 


131 Note differences in some M and Р values compared to EHTC V. We 


have corrected minor tabulation errors from that work, and have used a slightly 
different time range for averaging the MAD a, = —0.94 simulation. 

132 The Eddington rate is defined as Merga = Lgaa/€gaac?, where 
Leaa = АтСМт,с/от is the Eddington luminosity and we adopt an efficiency 
factor єваа = 0.1. Note that this assumed efficiency factor egag is distinct from 
the reported radiative efficiency є = L/Mc? measured from the simulations. 
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Figure 13. Average mass accretion rate (left panel) and jet power (right panel) for viable GRMHD models of M87” identified by selecting on total intensity data and 
jet power (blue, EHTC V), and when including polarimetric constraints (red). We estimate a mass accretion rate of M = (3—20) x 107^ Mo yr”', resulting in a 


radiative efficiency є < 1% (see EHTC V). The jet powers produced by our models are —-10*?? 


accretion rate is better constrained when including polarimetric information. 


(Di Matteo et al. 2003; see also Russell et al. 2015), and higher 
than that found from hybrid disk+jet models of the 
M87” spectral energy distribution (SED; Prieto et al. 2016). 
Our inferred jet powers of <10” erg s”! are toward the lower 
end of the observed range. In particular, the jet power measured 
at the location of Hubble Space Telescope (HST)-1 is 
10474 erg s^! (Stawarz et al. 2006), and Low-Frequency 
Array (LOFAR) observations suggest that а jet power of 104 
erg s ! was necessary within the last million years to inflate 
the observed radio lobes on scales of ^80 kpc (de Gasperin 
et al. 2012). 

Measurements of the accretion rate and the radiative efficiency 
can begin to constrain the microphysical plasma processes that 
heat electrons in M87", for example by inferring the fraction of 
the dissipated energy in the system that heats electrons, б. In 
axisymmetric, self-similar, hot accretion flow models, a system 
with M ~ 1075 Mgaa and a radiative efficiency є < 1% has а 
value of 6. in the range 0.1-0.5 (see Figure2 of Yuan & 
Narayan 2014). This range is consistent with that produced by 
simulations of turbulence and reconnection in the 3 ~ 1 regime 
(e.g., Rowan et al. 2017; Werner et al. 2018; Kawazura et al. 
2019). Future studies using simulations with self-consistent 
electron heating and radiative cooling (Section 6.3) can better 
constrain ó, and its dependence on local plasma parameters 
throughout the accretion flow and jet-launching region. 

We have assumed that all effects responsible for the 
appearance of the EHT polarized image of M87" are captured 
within the relatively small GRMHD simulation spatial domain, 
<102> rg. Goddi et al. (2021) developed a two-component 
model for the ALMA and image-integrated EHT data where 
each component is Faraday rotated by a different screen. The 
model demonstrates that the rotation measure of the compact 
component is unconstrained by the ALMA measurements 
alone, as the ALMA measurements are also sensitive to the 
Faraday rotation properties of the larger-scale component. In 
addition, the observed time variability in ALMA data (e.g., the 
RM sign change) can be explained by the observed EVPA 
variation of the compact core seen by the EHT. To produce the 
observed variability requires an RM of ~—6 x 10? rad m°. 
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—10* erg s^, and the jet efficiencies are ~5%-80%. The mass 


The ALMA data do not constrain the location or nature of this 
Faraday screen, except that it must be relatively close to the 
compact core, r < 10" r,. 

For our favored plasma parameters for M87", we expect 
substantial Faraday RM internal to the emission region itself, 
То, 2, 27, consistent with that measured from viable GRMHD 
images. In a model of uniform, external Faraday rotation this 
Faraday depth at 230 GHz would correspond to an RM of 
<10% rad m 2. Figure 14 shows that the apparent RMs 
measured from our GRMHD images span a wide range, 
often comparable to or larger than that inferred from the Goddi 
et al. (2021) two-component model (<10° rad m ?) For the 
low inclination angle of M87*, the apparent RM measured from 
GRMHD images is not a good tracer of the mass accretion rate 
(Moscibrodzka et al. 2017), and originates close to the 
emission region and well within the simulation domain (Ricarte 
et al. 2020, and Appendix B). The RM inferred from low- 
inclination GRMHD models of M87* can also vary rapidly and 
change signs (Ricarte et al. 2020), as seen in the ALMA-only 
data. As a result, the RM inferred from the two-component 
model in Goddi et al. (2021) is apparently consistent with the 
intrinsic properties of the GRMHD models studied here, 
without invoking an additional, external Faraday screen. At the 
same time, we cannot rule out that such an external screen 
could be present. Future EHT observations with wider 
frequency spacing can directly measure the resolved RM of 
the core and address this uncertainty. 


6.2. Electron Acceleration 


Magnetic reconnection, magnetohydrodynamic (MHD) tur- 
bulence and collective plasma modes in collisionless hot 
accretion flows likely result in nonthermal particle acceleration. 
While pioneering attempts have been made (e.g., Ball et al. 
2016; Chael et al. 2017; Davelaar et al. 2019), it is not yet 
known how to properly incorporate electron acceleration in 
global GRMHD simulations of hot accretion flows. 

We adopt an empirical approach to investigate the impact of 
nonthermal (accelerated) electrons on 230 GHz polarimetric 
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intrinsic polarization. 


images of M87'to quantify whether and how neglecting 
particle acceleration in our models affects our conclusions. We 
use a specific, but widely explored (e.g., Ozel et al. 2000; 
Markoff et al. 2001; Yuan et al. 2003), description for electron 
acceleration, namely that accelerated electrons add a power-law 
tail to the thermal distribution function. The power-law tail is 
described by the fraction of the thermal energy density in the 
power-law tail, 7, the power-law slope, p, and the maximum 
Lorentz factor of the accelerated electrons, 7y,,,,. The minimum 
Lorentz factor, yi, is calculated self-consistently by assuming 
a continuous transition between the thermal and power-law 
distribution functions (e.g., Yuan et al. 2003). In this model, the 
parameters p, т}, and 7, „, are constant across the accretion flow. 
We assume that y,,, — oo and we explore values of p = 3.5, 
4.5 and 7 = 0.01, 0.1. 

In Figure 15 we present linear polarization maps from two 
MAD models and one SANE model comparing purely thermal 
and hybrid electron distributions. Using a hybrid distribution 
function does not affect the structure of the EVPA map (62 
phase), but it changes the image-integrated and resolved linear 
polarization fractions. For example, in the MAD a, = —0.5 
(MAD а, = 0.5) model, with the selected hybrid parameters, 
the ||, Vnet and (|m|) ranges are 4.3%-4.6% (2.5%-3.8%), 
0.2596—0.3796 ((—0.5) to (—0.12)%), and 10.696—11.596 (12%- 
14%), respectively. Slightly larger deviations from the thermal 
model are measured in the SANE a, = —0.94 scenario, where 
the |т|һеь Vnet and (|m|) ranges are 2.296—4.196, —0.004% to 
0.31%, and 1496—20*6, respectively. 

However, fixing the accretion rate to that used in the thermal 
model results in an increased total intensity flux density when 
we add high-energy nonthermal electrons to our models. If 
instead we compare the models at fixed flux density, we need to 
reduce the mass accretion rate of the hybrid model. Therefore, 
generalizing the distribution function introduces order unity 
uncertainties in the inferred mass accretion rate, radiative 


IMInet 


Figure 14. Absolute value of RM vs. net linear polarization |m|ne for a subset of our EHT GRMHD library models explored in more detail in Ricarte et al. (2020). 
Closed symbols represent positive RM while open symbols represent negative RM, revealing significant time variability across the 2500 r,/c spanned by these 
snapshots. In gray, we plot our allowed region of |m|ne and bracket the range of core RM inferred from contemporaneous ALMA-only observations, 
2-100 x 10* rad m ? (Goddi et al. 2021). The dashed horizontal line demarcates the RM at which an EVPA rotation by т radians would have been observed between 
the 212 and 230 GHz frequency range used in the ALMA-only measurements, 1.05 x 10’ rad m ?. Despite large Faraday depths, a large fraction of these snapshots 
exhibit RMs consistent with simultaneous ALMA-only constraints. RM and [m|s« are anti-correlated, as larger Faraday depths lead to greater scrambling of the 
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efficiency, and jet power. The changes in the polarimetric 
observables in a given snapshot are also larger at fixed flux 
density. For example, in the MAD a, = —0.5 model, the |m|net 
increases from 4.796 to 6% when adding nonthermal electrons. 
As a result, in principle, polarimetric observables constrain the 
distribution function as well. 

Such constraints presumably depend on the details of the 
assumed particle acceleration scenario. Viable scenarios 
include hybrid electron distribution functions, or models where 
particle acceleration is a function of local gas conditions or 
magnetic energy density rather than fixed throughout the flow 
(e.g., Ball et al. 2016; Davelaar et al. 2019). More realistic 
particle acceleration scenarios could be considered using 
resistive GRMHD simulations (e.g., Ripperda et al. 2020). 


6.3. Coherently Polarized Forward Jet Emission 


As discussed above, some SANE retrograde model images in 
the library show coherently polarized features even when the 
Faraday depth through the entire emission region is large. The 
observed polarized flux in those cases originates on the near 
side of the midplane and is not scrambled from Faraday 
rotation along the line of sight. A similar effect might be 
possible if nonthermal electrons could be accelerated efficiently 
in the low-density, strongly magnetized funnel region in front 
of the black hole. 

It is beyond the scope of this Letter to evaluate whether or 
how such a model might be realized physically, e.g., whether 
any process could fill the funnel with high-energy electrons 
efficiently enough to produce the observed 230 GHz luminosity 
from the funnel alone. Instead we carry out one sample 
calculation of polarized emission from the funnel of a prograde 
d, = 0.94 SANE library snapshot. We assign a nonthermal 
energy density Unh = AUmag Wherever the magnetization с > 1, 
with a = 0.02 the fraction of the magnetic energy density Umag 
that is put into nonthermal particles. We calculate synchrotron 
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Figure 15. Sample polarization maps with varying electron distribution function. Columns display single snapshots from three selected models. Row 1 shows images 
with a thermal electron distribution function, as assumed in the standard EHT image library. Rows 2 through 5 are the same models but with emission from a hybrid 
distribution of electrons. Row 6 shows a hybrid model but the mass accretion rate of the model is adjusted to reproduce the same total intensity flux as the purely 
thermal snapshot. All maps are blurred with a 20 pas circular Gaussian. In all images, i= 17 deg (for negative a, models) or i= 163 deg (for positive a, model), 
with the black hole spin vector pointing to the left and away from the observer. 
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Figure 16. Sample library SANE a, = 0.94 Rnigh = 160 snapshot (left panel) and with power-law emission from nonthermal electrons added in c > 1 regions (in the 
funnel near the pole, middle panel) and in the funnel only outside of a radius ra; = 6 r, (right panel). The library model is heavily depolarized due to Faraday rotation. 
Nonthermal radiation from the forward jet is coherently polarized; these images look qualitatively similar to polarimetric images of the forward-jet dominated semi- 
analytic models of Broderick & Loeb (2009). Even a small total flux contribution can increase the net and image-averaged linear polarization fractions to lie within the 
observed range. However, in this example ће БУРА patterns (8 phase) of the TH--PL images remain inconsistent with M87” data. 


radiation from a pure power-law distribution of electrons with 
Ymin = 100 and p = 3. 

Figure 16 compares the original thermal snapshot with two 
realizations of this hybrid thermal4-nonthermal funnel emission 
models. In the purely thermal case, Faraday rotation depo- 
larizes the emission at the EHT beam scale, producing low 
fractional polarization across the image that is inconsistent with 
EHT observations of M87*. Adding power-law electrons in the 
funnel produces coherent linearly polarized emission. When we 
assume Unth X Umag (middle panel), the power-law emission is 
concentrated close to the black hole and lensed into a ring 
(Dexter et al. 2012). The weak forward jet component is 
strongly polarized but lies inside the observed ring, and is thus 
potentially inconsistent with the EHT total intensity and 
polarimetric image. In the right panel, we exclude nonthermal 
emission from inside a radius ra = Org. Both nonthermal maps 
are consistent with our cuts on net and image-average linear 
polarization fraction, |т| and (|m|) However, both are 
inconsistent with the observed EVPA pattern of M87" (i.e., the 
б» phase). 

For this example, we assume a plasma of protons and 
electrons rather than e*/e' pairs. The latter are presumably 
more likely to form in the funnel (Moscibrodzka et al. 2011; 
Chen et al. 2018; Levinson & Cerutti 2018; Anantua et al. 
2020; Crinquand et al. 2020; Wong et al. 2021), and have 
different circular polarization properties. Future observations 
that constrain the resolved circular polarization structure could 
potentially discriminate between pair and electron-ion plasmas 
in the emitting region. At longer wavelengths and larger scales, 
the limb-brightened jet structure of M87 (e.g., Walker et al. 
2018) also suggests that the radiating electrons are not 
concentrated inside the funnel as modeled here. 


6.4. Radiative Models 


Our GRMHD images use the parameterization of 
Mo$cibrodzka et al. (2016) to model the electron and ion 
temperatures given the total gas temperature from a simulation. 
In this prescription, the electron-to-ion temperature ratio is a 
function entirely of the local plasma 5. This functional form 
(Equation (15)) captures the general behavior seen in many 
simulations of electron heating in turbulent or reconnecting 
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collisonless plasmas; namely, the electron heating is more 
efficient (and thus the temperature ratio is closer to unity) when 
B « 1 (e.g., Howes 2010; Rowan et al. 2017). However, the 
actual distribution of 7T, in a hot accretion flow reflects the 
balance of heating, cooling, and advection of hot electrons 
throughout the system. Furthermore, the GRMHD simulations 
in the library considered here do not include radiative cooling. 
Our passing models for M87* favor a radiative efficiency of 
€ ^» 1% (Section 6.1), however, and we may begin to worry if 
cooling is dynamically important in M87”. 

To assess these uncertainties, it will be useful to compare the 
results in this work with results from simulations performed 
with radiative GRMHD codes. These codes typically use either 
the MI closure method (e.g., Sadowski et al. 2013; McKinney 
et al. 2014; Sadowski et al. 2017) or a Monte Carlo approach 
(e.g., Ryan et al. 2015) to track radiation and its interactions 
with the plasma near the black hole. In addition to the effects of 
cooling on the gas temperature, these codes can also evolve the 
separate electron and ion temperatures under the influence of 
cooling and different subgrid prescriptions for the electron 
heating efficiency (e.g Ressler et al. 2015, 2017; Chael et al. 
2018, 2019; Ryan et al. 2018; Dexter et al. 2020). 

In Figure 17, we show a comparison of the temperature ratio 
T;/T, obtained directly from the two MAD radiative simula- 
tions of M87* in Chael et al. (2019; right column), with the 
temperature ratio obtained from the same simulations using 
Equation (15) with Riow = 1, Rpigh = 20 (left column). The two 
rows show simulations using different underlying models for 
electron heating: the top row (H10) uses the turbulent heating 
prescription of Howes (2010), and the bottom row (R17) uses 
the reconnection heating prescription of Rowan et al. (2017). 
The simulation data in both rows is averaged in time and 
azimuth. 

Figure 17 shows that the temperature ratios obtained from 
the electron-ion evolution and from the postprocessing 
prescription both transition from smaller to larger values when 
moving from the jet/funnel region (low 8) to the disk (high 5). 
In simulation H10 (top row), T;/T, & 1 in the funnel, which 
well matches the result from Equation (15) with Riow = 1; in 
simulation R17, electrons are cooler in the funnel (T;/T, ~ 5), 
and even colder in the disk-jet interface directly outside the 
о = 1 contour (T;/T, e 15). The cooler electrons in the funnel 
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Figure 17. Comparison of the temperature ratio obtained directly from two 
radiative GRMHD simulations with the values obtained from the postproces- 
sing model used in this work. The top and bottom rows show, respectively, 
time- and azimuth-averaged data from the two radiative MAD simulations of 
M87* presented in Chael et al. (2019). The top row shows data from simulation 
H10, where electrons are heated by the turbulent heating prescription of Howes 
(2010), and the bottom row shows the R17 simulation, where electrons are 
heated by the magnetic reconnection prescription of Rowan et al. (2017). The 
left column shows the spatial distribution of the ion-to-electron temperature 
ratio obtained from these simulations by applying our /j-dependent post- 
processing equation to the total gas temperature (Equation (15) with Rio, = 1, 
Rhigh = 20). The right column shows the temperature ratio obtained directly 
from the independently evolved electron and ion species in the radiative 
simulations. The white contour shows the c = 1 surface. 
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regions of R17 reflect the decreased efficiency of low-@ 
electron heating in the Rowan et al. (2019) model than in 
Howes (2010). In Section 5.4, we note that the results of joint 
scoring (but not simultaneous scoring) favor Riow = 10 models 
over Riow = 1 models. Model R17 shows that some radiative 
models can naturally produce Т,/Т, > 1 in low-@ regions. 
However, the funnel value of T;/T, ~ 5 in this simulation is in 
between the two cases Riow = 1, Riow = 10 that we considered 
in the image library. 

While the disk electrons are cooler than the jet electrons in 
both radiative simulations shown in Figure 17, T;/T, in the disk 
of model H10 increases at small radii. In contrast, model R17 
better matches the postprocessing model prediction of a 
roughly constant disk temperature ratio. (Note that the 
230 GHz emission is entirely produced at radii r  10r,.) 

In addition to the MAD simulations from Chael et al. (2019) 
in Figure 17, we have also checked the temperature ratios in the 
radiative SANE simulations of Ryan et al. (2018). In these 
simulations, the average temperature ratio in the EHT 230 GHz 
emission region can also be roughly approximated by the 
Mo$cibrodzka et al. (2016) prescription, with А = 1 and 
Ryign in the range Rhign ~ 10—20. 

Our preliminary results indicate that while some important 
features of the temperature-ratio distributions produced in 
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radiative simulations can be described by the Riow, Rnign model 
(Equation (15), the current postprocessing model cannot 
capture all of the behavior produced in radiative simulations. 
A more detailed comparison is left for future work. 


7. Predictions 


We have identified a subset of a large parameter space of 
GRMHD models that is consistent with constraints derived 
from current EHT total intensity and polarimetric observations 
of M87". The models that pass our constraints on the 
polarimetric structure and jet power from М87*аге all 
magnetically arrested (MAD) accretion flows. Here we make 
predictions for testing our interpretation with future 
observations. 


7.1. Repeated Observations 


Repeated EHT observations of M87* at 230 GHz (both in 
total intensity, e.g., Wielgus et al. 2020, and in linear 
polarization) will continue to constrain the model parameter 
space. Figure 18 shows the time evolution of 3 amplitude and 
phase for 200 snapshots of three viable library models: MAD 
ax = —0.5, Row = 10, Rhigh = 20; MAD аж = +0.5, Riow = 10, 
Rhign = 80; and MAD ay = +0.94, Riow = 10, Rhigh = 80. The 
observer inclination was 17 deg and 163 deg for the retrograde 
and prograde models, respectively. 

Both quantities show variations on timescales from days to 
months. The phase and amplitude of 3, should change over the 
course of a week of observations. In EHTC VII, we observe 
changes in the the (5> amplitude and phase over the week of 
observations in 2017, and use the results from two epochs to 
define our acceptable parameter ranges. Figure 18 suggests that 
occasionally the observed changes in (5 on ~week timescales 
can be much more dramatic than we observe in 2017, with 
variations in (5 phase of 90 deg for some models on short 
timescales. 

The scatter in both quantities on longer “month timescales is 
much larger than the uncertainty range derived from the EHT 
2017 measurements. If our passing GRMHD models accurately 
describe the 230 GHz emitting region in M87', future EHT 
observations should detect variability in the polarization 
structure. According to current models, the time-averaged 62 
amplitude and (|) should remain similar to the current values 
for prograde spin models, and tend toward larger values for 
retrograde spin models. For high-prograde spin (or many 
SANE models), the 8> phase should on average be closer to 
zero than we observe in 2017. 


7.2. Future Observations at 260 and 345 GHz 


In selecting models, we have focused on metrics corresp- 
onding to salient features of the data. We have not attempted to 
compare models in detail to specific features of the recon- 
structed polarimetric images, most notably the apparently 
depolarized bright patch in the eastern part of the image 
(Figure 1). We do note that such depolarized features occur in 
many of our library images, particularly in MAD models with 
Riow = 10. If the eastern patch in the 2017 image is depolarized 
due to Faraday rotation, it may be possible to tell with future 
higher-frequency observations. Figure 19 shows images of two 
of the favored MAD models at the current observing frequency 
of 230 GHz and at two additional frequencies planned for the 
future EHT observing campaigns. In addition to internal 
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Figure 18. Amplitude (left panel) and phase (right panel) of £5 as a function of time for three viable GRMHD library models identified here (points, all with 
Коу = 10) compared to ranges measured from ЕНТ 2017 M87* data (gray shaded region). The dashed lines show the median values for each model. The retrograde 
spin model predicts higher (2 amplitude in future observations. In the high-prograde spin model, the median (5 phase is closer to zero than the observed range in 2017. 
Changes in both quantities occur on timescales of weeks to months, and should be apparent in future EHT data sets. 
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Figure 19. Random snapshot of MAD models with а, = —0.5 (top row) and with a, = 0.5 (bottom row) at the current EHT frequency of 230 GHz and two higher 
frequencies of 260 and 345 GHz, which are planned for the future ЕНТ observations. All images are convolved with a 20 pas Gaussian. In all images the black hole 
spin vector is pointing to the left and away from the observer. In all cases, the ring fractional polarization increases slightly with frequency. The EVPA pattern, as 


measured by the (3, is similar at all three frequencies. 


Faraday rotation, the sense of the EVPA pattern may also be 
subject to a net, coherent rotation due to external Faraday 
rotation. At higher frequency, Faraday rotation is suppressed 
and EHT observations will see the intrinsic magnetic field 
pattern more clearly. 

For a snapshot from the MAD a, = —0.5 (Rhigh = 160 and 
Riow = 1) model, Figure 19 shows that the |m|,4 and (|m|) 
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values are predicted to increase with frequency. The 230, 260, 
and 345 GHz net EVPAs are —77, —70, and —82 deg, 
respectively, corresponding to an (apparent) rotation measure 
RM ~ 1 x 10 rad m ? between 230 and 345 GHz. The net 
circular polarization |v|net remains small and nearly constant 
with frequency; it is 0.42%, 0.35%, 0.32% for 230, 260, and 
345 GHz, respectively. A similar trend is observed in a MAD 
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dy = 0.5 (Rnigh = 80 and А, = 10) model. The image |m|net 
and average polarization (|m|) are again expected to increase 
with frequency. The corresponding net EVPAs are —38, —40, 
and —30 deg, corresponding to an apparent rotation measure of 
RM ~ —1 x 10? rad m ?. The net circular polarization fraction 
[уос remains roughly constant and close to zero, 0.33, 0.06, 
and 0.2% from low to high frequency. 

Both of these models display similar EVPA structure at all 
three frequencies, indicating that in this example the net EVPA 
pattern is due to magnetic field structure rather than coherent 
Faraday rotation. Future multi-frequency observations will be 
able to infer the core RM and intrinsic EVPA pattern set by the 
near-horizon magnetic fields. 


8. Conclusions 


The EHT has produced resolved polarized intensity maps in 
the near-horizon region around the supermassive black hole in 
M87. Taken together with image-integrated data from simulta- 
neous observations with ALMA, these images constrain the 
space of accretion flow and jet models used to interpret the 
EHT total intensity image with broad implications for jet 
launching near a black hole event horizon. Here we summarize 
the main results of that analysis. 


1. We interpret the depolarization seen in EHT images as 
the result of beam depolarization due to Faraday rotation 
internal to the emission region. In the context of one-zone 
models and combined with the size and brightness 
temperature of the total intensity image, we estimate an 
average emission region plasma density of ne~ 
1047 cm ?, magnetic field strength of B ~ (1-30) G, and 
T, = (1-12) x 10'? K. 

2. The net БУРА pattern of the M87” polarization maps is 
predominantly azimuthal. In the context of semi-analytic 
models with imposed, idealized magnetic field geometry, 
such a pattern can be reproduced using a significant 
component of poloidal (radial and/or vertical) magnetic 
field. The presence of such magnetic fields in a rotating 
fluid would imply that the magnetic fields are dynami- 
cally important. However, significant Faraday rotation 
may be present and it is not clear whether the observed 
EVPA pattern can be interpreted in terms of magnetic 
field structure alone. 

3. To capture the effects of realistic magnetic field structure, 
plasma conditions, and Faraday rotation and conversion, 
we have compared salient observables to a large library of 
images from the GRMHD simulation library of EHTC V. 
The observables are the net circular polarization fraction 
constrained by ALMA (|v|nc), the net and image- 
averaged linear polarization fraction measured by the 
ЕНТ (| nee and (|m|)), and the m = 2 coefficient of a 
Fourier expansion of the azimuthal EVPA pattern (62). 
Of these, 9; is the most constraining metric. 

4. The model-scoring procedures disfavor most models 
from the GRMHD image library from polarimetric 
observations alone. Many weakly magnetized (SANE) 
models are too depolarized, or show an EVPA pattern 
that is too radial. Many strongly magnetized (MAD) 
models are too coherently polarized. The polarization 
fraction is generally set by the Faraday rotation depth 
close to the emission region. MAD models more 
frequently produce azimuthal EVPA patterns, as expected 
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for magnetic field structures that include a significant 
poloidal field component. Combined with a conservative 
lower limit on the jet power of M87, only strongly 
magnetized (MAD) models remain viable. We use those 
remaining models to estimate the mass accretion rate 
onto the central supermassive black hole as M — 
(3-20) x 107^ Mo yr-!. The average plasma parameters 
found from GRMHD images are in good agreement with 
those inferred from one-zone models. 

5. The model space considered in this Letter is incomplete, 
and systematic uncertainties remain a challenge. While 
the radiative efficiency that we find is relatively high, we 
consider only non-radiative GRMHD models. We do not 
consider GRMHD models with misalignment between 
the disk and the black hole angular momentum. We also 
only consider one parameterization for determining the 
electron distribution function from the simulation data. Of 
these three major areas of uncertainty, we have explored a 
small sample of alternative models for determining the 
electron distribution function, including both alternative 
prescriptions for electron heating in strongly magnetized 
regions and including a nonthermal component. The 
quantitative estimates of mass accretion rate and jet 
power found here depend on the assumed electron 
distribution function and are uncertain at the order unity 
level. The alternative electron distribution functions 
considered here do not change the main finding that 
MAD models with dynamically important near-horizon 
magnetic fields appear more viable for explaining the first 
polarimetric EHT observations of M87”. 

6. Our favored models show time variability in the 
polarization metrics used here. The median values found 
at several epochs should be sufficiently well measured to 
distinguish between the current retrograde and prograde 
spin models. At higher frequencies of 260 and 345 GHz, 
weaker Faraday effects should result in an increased 
degree of polarization. Measurements of the EVPA 
pattern at higher frequencies can distinguish between 
Faraday rotation along the line of sight and the imprint of 
the underlying magnetic field structure. Continued 
imaging with the EHT and advances in radiative and 
nonthermal theoretical models will further constrain the 
electron distribution and magnetic field structure in the 
jetlaunching region near the supermassive black hole 
event horizon in M87. 
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Appendix A 
Relationship between 62 Coefficient and E- and B-modes 


The (5 coefficients of the azimuthal decomposition of the 
complex linear polarization P = О + iL used in Section 5 are 
directly related to the decomposition of polarization fields into 
E- and B-modes familiar from cosmology. In this Appendix, 
we illustrate that relationship and demonstrate that the 
information in the image-space decomposition of GRMHD 
library images in Section 5 can also be accessed directly in 
calibrated visibility domain data sampled on EHT 2017 base- 
lines, provided the data are accurately phase calibrated. 
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A.1. Definitions 


In defining flat sky E- and B-modes, we follow the 
conventions of Kamionkowski & Kovetz (2016), Section 4.1 
(up to factors of 4/2). E- and B-modes are naturally defined in 
the visibility space sampled by an interferometer. For a baseline 
vector u with a magnitude и and PA Ө, the E and B visibilities 
are related to the Stokes visibilities 9 and M by a rotation of 
20 in the Fourier plane 


E(u, 0) _ | cos 20 m Au, 0) 

Bu, 0) —sin20 cos20 | tu, D) 
In real space, the E- and B-mode fields are analogous to the 
gradient and curl of the polarization tensor: 


(Al) 


VE = д,0Р,, V?B = єадьд.Р,ь, (A2) 


where €ac is the 2D Levi-Cevita symbol and the polarimetric 


tensor is 
Q u 
Pa = : 
? E A] 


P transforms as a trace-free tensor under rotations; for a 
rotation matrix R(o) that rotates the coordinate axes by an 
angle o, P — R(o)PR'(o) (equivalently, the complex field 
P — Pe”), While the values of the Q and U images depend 
on the choices of coordinate axes, the real space E- and B-mode 
images are coordinate-independent scalars. 


(A3) 


A.2. Relationship between (E, B) and 3, Coefficients 


Consider a linearly polarized image P = Q + iU in 2D 
image-domain polar coordinates (p, Фф). We can expand the 
image in a multipole series: 


P(p,d) = № У` Bf, (oem. 


m=— оо 


(A4) 


In the decomposition of Equation (A4), Jo is the total flux 
density of the Stokes Z image, the 8, coefficients are complex, 
and the radial envelope function f,,(p) is normalized so that 


max 


(A5) 


р, 
2т f 
Pmin 


Jn (p)p ар = 1. 


The Øm coefficients defined in this way then correspond to 
those defined in Equation (9): 


] рх 
Bm = i f 


Pmin 


f P(p, Ф)еіп®р dp dé. (A6) 


In particular, бо is the image-integrated complex fractional 
polarization, and £5 encodes the same information on the 
gradient and curl of the polarization field that is available in the 
E- and B-modes. Note that because P is a complex image, in 


general Bm = B* m 
The Fourier transform of P(p, ф) is 


oo 


P(u, 0) = 2710 PA i™ bme"? F (U), 


m=— oo 


(A7) 
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Figure 20. Schematic guide to EVPA patterns for ring images with power only іп the 3, mode in an azimuthal decomposition, and the corresponding signs of the Ё 


and B mode visibilities. 
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Figure 21. Distributions of the amplitude (left panel) and phase (right panel, in degrees) of the complex quantity ев = — Єв — iBg computed from simulated data on 
ЕНТ 2017 baselines for the full GRMHD image library considered in this work. The distributions are broadly consistent with the (5 results measured in the image 
domain in Figure 9, illustrating the relationship between the 3, metric and average Ё and B mode visibilities. 


where F,(u) is the mth order Hankel transform of the radial 
function f,,(p): 


Pmax 
F0) = [^ „(о).Һ(2три)р dp. (A8) 
Prin 


The Fourier transforms of Q and UY (the linear Stokes 
visibilities) are then 


00 


Au, 0) = 21lo У) i" Re[B,e"^F, (u)], 
(и, 0) = 2710 > i7" Im[8,,e"*F, (u)], (A9) 


and from the definition in Equation (A1) the E- and B-mode 
visibilites are 


оо 


Eu, 0) = 2719 Уу i" Re [Bme 99 E, (u)], 
B(u, 0) = 271) У) i""Im[B,e"-29?F,(u). (А10) 


From Equation (А 10), ме сап see immediately that ап image 
with real (5 and all other 8,22 — 0 is a рше E-mode; for 
instance, if 8 = 1 (radial polarization vectors), Ё < 0. If 
9; = —1 (toroidal polarization vectors), Ё > 0. Similarly, if an 
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image has purely imaginary ® (and all other 8,2 = 0) it is a 
рше B-mode; if e.g., 3, = i (right-handed helical polarization 
vectors), then B < 0, and if = —i (left-handed helix), 
B > 0. Figure 20 illustrates the values of the (5> and the 
corresponding signs of ће Ё and В visibilities for these four 
azimuthally symmetric cases. 


A.3. É and B Distributions of GRMHD Library Images 
Starting with Equation (A10), and integrating over the 
baseline angle 0, only the 3, mode survives. 


J. i E(u, 0)d0 = —2rh Re [82 Е›(и)] (A11) 


2n 
f B(u, 0)d0 = —2rh Im [85 (W1. (A12) 

Thus, if we have calibrated visibility measurements of 0 and 
U (and thus Ё апа B), we can measure the (5> mode by 
averaging É and B in visibility space. 

To illustrate this connection between E- and B-mode 
visibilities and the 25 coefficient, we compute the following 
averages on images in the GRMHD library 


ne (Re[É])q v») u ( Re[B])q, v») 
R= es R= y 


= = (А13) 
(а) (А Na) 
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where we take the average only over (u, v) points sampled by 
the ЕНТ in 2017, including conjugate baselines. As atmo- 
spheric phase errors and D-term miscalibration will affect 
the phase of the Ё and B visibility and thus the results for the 
averages, we perform this test on synthetic data from the image 
library using perfectly calibrated data with no noise. 

Because we include conjugates in the average over all data 
points, the average of the imaginary part is zero. We perform 
the averaging only over а и, v range (1, 10]СА in order to 
remove any effects from large-scale structure, which may have 
a different net sense of polarization than the resolved 
emission ring. 

We finally combine £g and Bg in a complex quantity: 


beg = —Ex — iBsg, 


where the negative signs are chosen to match the angle 
convention for 6 in Equation (9). In Figure 21 we show 
histograms of the magnitude and angle of the (eg for 
comparison with the 5 histograms in Figure 9. The 
distributions broken down by model type (MAD/SANE, 
prograde/retrograde) reproduce the general behavior of (5 
amplitude and phase from the image-domain calculations in 
Figure 9, although the normalization of the gg amplitude is 
different. Comparing Figure 9 with Figure 21, it is apparent that 
all of the essential information on the EVPA structure used in 
this Letter can, in principle, be extracted from EHT visibilities 
without image reconstruction. However, because phase and 
amplitude calibration of the EHT visibilities is necessary for 
extracting the E- and B-modes from the visibilities, modeling 
the source structure in the image domain would remain a 
necessary part of the analysis even if we were to use £g and Bg 
instead of 5. 


(А14) 


Appendix В 
Faraday Rotation in GRMHD Models of М87* 


As linear polarization travels through magnetized plasma, 
Faraday rotation shifts its EVPA by 7,, /2 radians. If 7, > 1, 
as is the case for most of our models (see Figure 7), Faraday 
rotation can, in principle, scramble otherwise observable 
polarimetric signals. In this section, we explore in more detail 
the sources of Faraday rotation in our models, and demonstrate 
that observable linear polarization signals can, in fact, exist in 
models with 7,, > 1. This is because Faraday rotation occurs 
co-spatially with the emission and should not be conceived of 
as a purely external screen. 

Ricarte et al. (2020) studied the resolved Faraday rotation 
properties of a subset of the same models used in this work. 
Figure 14 shows their inferred |RM| versus |m|;« for those 
images. For each model, 11 snapshots spaced between 7500 
and 10000 r,/c are included. Each of these models was found 
to pass the constraints of EHTC V. Snapshots with positive 
КМ5 are plotted with filled symbols, while those with negative 
RMs are plotted with open symbols. In gray, we overplot the 
allowed range of |m|net as well as the range of RM for the core 
region inferred from simultaneous ALMA-only observations. 

Despite the large Faraday depths of these models, many of 
them are capable of producing RMs that are consistent with the 
observed data. RM and |m|net are anti-correlated, as expected, 
because a greater amount of Faraday rotation should both 
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increase the RM and cause a greater amount of scrambling of 
the polarized emission. Note that the RM varies by orders of 
magnitude and even flips sign over time in these models. This 
is due to the summation of time-variable regions with 
significantly different and even oppositely signed Faraday 
rotation depths, which also contributes to highly non-AX 
evolution of the EVPA with wavelength. 

One potential source of uncertainty is the limited volume of 
our GRMHD simulations. In our image library, the radiative 
transfer equation is solved within a radius of only 50 or 100 r, 
depending on the model, while in principle significant Faraday 
rotation may occur at much larger radius. Figure 22 
demonstrates that for the gas distributions studied, most of 
the Faraday rotation occurs at radii much smaller than the outer 
domain. Five example models are visualized here: (a) MAD 
ax = +0.94 Къв = 20, (b) SANE аж = +0.5 Rhigh = 1, (c) 
MAD a, = —0.5 Ryign = 160, (d) SANE а, = +0.5 Кыһ = 
160, and (e) SANE a, = 0 Аһ = 80. The brightness of each 
pixel scales with the total intensity (intentionally saturating 
0.396 of the pixels), while the color of each pixel scales with 
the intensity-weighted Faraday depth 7,,;, as shown in the 
colorbar. ту is distinct from т, because it is intensity 
weighted along each ray, such that in each pixel 

Tt = F floras. (В1) 

As also shown in Figure 7, these models span а wide range 
of Faraday depths. Typically, SANE models and models with a 
larger Rpign have larger Faraday depths than MAD models and 
those with smaller Къы. SANE models require a larger 
accretion rate to reproduce the total intensity of M87", which 
increases the amount of Faraday rotating material. Meanwhile, 
increasing Ay, lowers the temperature of the midplane by 
construction. This makes Faraday rotation more efficient, and 
also requires a larger accretion rate to compensate for the lower 
electron temperatures (see Mo$cibrodzka et al. 2017, for ап 
extended discussion). 

In Figure 23, we confirm that the linear polarization parameters 
used in this study are not strongly evolving at the outer simulation 
domain. Here, we provide polarization maps and the linear 
polarization parameters for these models blurred with a Gaussian 
beam with a FWHM of 20 pas. The rows display models with 
outer integration radii of 10, 20, 40, and 50 r,. For the three left- 
most models, there is little difference between images constructed 
with max = 10 r, and those with Max = 50 r,, echoing our 
previous findings that there are small fractional differences in the 
Faraday depth between these scales. Faraday rotation thick 
models (d) and (e) show modest differences between images 
calculated with outer boundaries of 10 r, and 50 r,. Those images 
also appear to be converged by 40 ry. 

Some of our models produce observable polarimetric 
signatures despite 7,, being large enough to potentially 
depolarize all of the emission. This apparent contradiction is 
resolved by the fact that not all emission is Faraday rotated by 
the same amount. Because Faraday rotation occurs co-spatially 
with the emission, instead of as an external screen, there can 
exist emission traveling on Faraday-thin paths to the camera 
even in models where 7,, > 1 when integrated along the entire 
geodesic. In Figure 24, we illustrate this phenomenon by 
splitting these images into the emission originating in front of 
or behind the midplane. Here, models are plotted as in 
Figure 22 with љах = 50 r,. Emission originating from the 
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[may = 20 rg 


y [uas] 


A 


Toi =4.39x 10-1 TAN 29807107  1,,7123.57X10'  15,25./06x10^ Topi = TSL x 105 


Imax = 50 rg 


y [uas] 


P d 


To,174.37 X10! T, М. 29109 т,,,=3.56х 101  15,,"25.82X10* To= 1:32 Xx 10? 
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Figure 22. Intensity-weighted Faraday depths visualized with five example models: (a) MAD a, = +0.94 Кмьһ = 20, (b) SANE a, = +0.5 Кмьһ = 1, (c) MAD 
ay = —0.5 Ryigy = 160, (d) SANE a, = +0.5 Аһ = 160, and (e) SANE a, = 0 Аһ = 80. The brightness of each pixel scales with the total intensity (intentionally 


saturating 0.3% of the pixels), while the color indicates the intensity-weighted Faraday depth, т, ;. In the top row, the maximum integration radius is set to 20 r,, 


while in the bottom row, the maximum integration radius is set to 50 r,. We find very little difference, confirming that our results should be insensitive to the outer 
radius of the simulation domain. 


back half of the simulation domain exhibits larger values of Emission from in front of the midplane that is observed within 
Тр, AS it must travel through more Faraday rotating material. the photon ring has a much larger T, than the rest of the 
Notice the clearly Faraday-thin regions in panels (c) and (e), image because those geodesics pass through the midplane and 
despite the enormous values of (7,,,7) for the image overall. around the black hole through Faraday-thick material. 
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Fractional Polarization m 
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Figure 23. Polarimetric images of five example models: (a) MAD a, = +0.94 Rhigh = 20, (b) SANE a, = 40.5 Rug, = 1, (c) MAD a, = —0.5 Rhign = 160, 
(d) SANE a, = 10.5 Ruign = 160, and (e) SANE a, = 0 Rnigh = 80. The maximum integration radius is set to 10, 20, 40, and 50 r, in each row from top to bottom. 
For models (a)-(c), there is little difference between images computed with 7,4, = 20 г, and those computed with max = 50 r. The significantly more Faraday-thick 
models (d) and (e) appear to be converged by a radius of 40 г,. 
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Figure 24. Faraday depth visualizations as in Figure 22, but with emission origin split into the front and back halves of the simulation domain. Fmax = 50 r for all of 


these images. Because the Faraday rotation occurs co-spatially with the emission, radiation originating 


from the front half of the simulation domain has smaller Toy 


than emission originating from the back half. In panels (c) and (e), notice the Faraday-thin (т, < 1) regions (purple) in the front half images even though 7,,,; > 1 


for the model overall. 


Appendix С 
Distributions of Theory Metrics for Each Model 


Figures 25-29 show distributions of the metrics used to 
score models, split out for each model individually. 

Prograde SANE models show rapidly decreasing |m|;« with 
increasing Rhign (Figure 25) and are significantly depolarized 
when Rhign 10. This behavior was previously demonstrated 
by Mo$cibrodzka et al. (2017). The accretion flow electron 
temperature decreases with increasing Ryign, increasing the 
strength of Faraday rotation while also concentrating the 
emission at high latitudes behind the black hole (see 
also EHTC V). The emission is then depolarized when 
traveling through the Faraday-thick midplane plasma. 

Retrograde SANE models, however, show nearly the 
opposite behavior, with depolarization maximized for 
Ryign = 1. At larger values of Ryigr, linearly polarized emission 
appears on the near side of the midplane, producing coherent 
linear polarization structure that is not Faraday depolarized. 

MAD models at all spins show a mild degree of 
depolarization with increasing Rpignh. The accretion flow 
electron temperature remains high even for large values of 
Rhign, as much of the plasma has 3 2 1. 

Similar qualitative behavior is seen in (|m|) (Figure 26) and 
the amplitude of (5 (Figure 27). However, those quantities 
show less time variability (narrower distributions) than is seen 
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in |m|net- AS a result, observed ranges of those values are more 
constraining. In particular, the MAD models show consistent 
offsets where ([m|) and || are lower for Ао = 10 than 
Riow = 1 models. Some spin dependence is also apparent, with 
high-prograde spin usually corresponding to the highest 
degrees of ordered polarization. 

When the (> amplitude is not strongly suppressed (e.g., by 
Faraday rotation), the (5 phase distributions are related to 
intrinsic magnetic field structure (e.g., Figure 3 and PWP). 
Prograde spin, Аһ = 1 SANE models and retrograde spin, 
large Rpign SANE models both show radial EVPA patterns, 
resulting in £5 phase distributions near zero. MAD models 
show spin-dependent ( phase distributions for low values of 
Rygs ranging from spiral patterns (ZG, = —90 дер) for 
retrograde spin to more radial patterns at high-prograde spin. 
The patterns are relatively constant functions of Rpigh and Riow, 
although with some shift of MAD prograde distributions to 
twistier EVPA patterns, particularly for Rigy = 10. 

Most models show distributions of ve; centered on zero, 
near the observed range (Figure 29). MAD models generally 
show low circular polarization fractions, while heavily 
depolarized SANE models (retrograde low Rhigh, prograde 
large Rpigh) tend to show larger IVnet| than observed, which сап 
be explained by stronger Faraday conversion in the emission 
region. 
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Figure 25. Distributions of net polarization fraction |m|net for all models. MAD and SANE simulations are shown in the left and right panels. Black hole spin а, varies 
along the x axis, Rhign varies in each row, and the distributions at R,,,, = 1 and 10 are shown in red and blue in each case. 
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Figure 26. Same as in Figure 25, but for the image-averaged polarization fraction (||). 
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Figure 27. Same as in Figure 25, but for ће amplitude of (3, |9]. 
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Figure 28. Same as in Figure 25, but for the phase of 62. 
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Figure 29. Same as in Figure 25, but for the net circular polarization fraction Vnet- 


Appendix D 
Detailed Model Scoring Results 


In Table 3 we provide a summary of the number of images 
for each model that fall within the observed range of each 
individual theory metric (used in the joint scoring procedure) 
and within the observed ranges of all metrics simultaneously 
(used in the simultaneous scoring). Boldfaced type is used for 
models that are deemed viable by one of the scoring systems. 
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For simultaneous scoring, a viable model contains at least one 
image that simultaneously satisfies all constraints. For joint 
scoring, a viable model has a joint likelihood >1% that of the 
maximum found across all models. We also provide a summary 
score—“pass” indicates a model that satisfies the polarimetric 
constraints according to either scoring procedure, as well as the 
jet-power cut of Piet > 102 erg s! (EHTC V). 
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Table 3 
Scoring Results for the Models Used Here 
Flux ax Riow Кв М.з Радо № № Мт) Nga Narg б Ми Summary 
SANE —0.94 1 1 3.79 1.19 7 578 0 0 126 0 fail 
SANE —0.94 1 10 16.30 5.11 504 95 0 0 0 0 fail 
SANE —0.94 1 20 19.04 6.00 482 24 0 0 0 0 fail 
SANE —0.94 1 40 22.88 7.26 460 27 0 0 0 0 fail 
SANE —0.94 1 80 29.18 9.19 375 31 0 0 0 0 fail 
SANE —0.94 1 160 39.32 12.40 301 41 0 0 0 0 fail 
SANE —0.94 10 1 6.24 1.96 0 506 0 0 185 0 fail 
SANE —0.94 10 10 48.18 15.16 263 26 28 12 0 0 fail 
SANE —0.94 10 20 55.54 17.47 256 32 27 6 0 0 fail 
SANE —0.94 10 40 66.19 20.82 256 36 34 7 0 0 fail 
SANE —0.94 10 80 83.97 26.42 251 30 39 11 0 0 fail 
SANE —0.94 10 160 115.94 36.47 232. 27 49 15 0 0 fail 
SANE —0.5 1 1 2.48 0.04 0 407 0 0 143 0 fail 
SANE —0.5 1 10 18.36 0.30 51 0 10 15 0 0 fail 
SANE —0.5 1 20 21.37 0.35 78 0 12 28 0 0 fail 
SANE —0.5 1 40 24.80 0.41 61 2 10 16 0 0 fail 
SANE —0.5 1 80 30.28 0.50 63 13 5 12 0 0 fail 
SANE —0.5 1 160 39.73 0.66 80 37 D 7 0 0 fail 
SANE —0.5 10 1 3.75 0.06 0 491 0 0 134 0 fail 
SANE —0.5 10 10 64.76 1.06 9 3 0 0 13 0 fail 
SANE —0.5 10 20 84.43 1.39 107 0 0 0 10 0 fail 
SANE —0.5 10 40 92.84 1.53 95 0 0 0 13 0 fail 
SANE —0.5 10 80 107.39 1.77 70 0 0 0 9 0 fail 
SANE —0.5 10 160 134.03 2.20 48 0 0 0 11 0 fail 
SANE 0.0 1 1 0.89 0.00 555 513 90 96 38 0 fail 
SANE 0.0 1 10 17.40 0.00 78 174 12 29 11 0 fail 
SANE 0.0 1 20 31.92 0.00 334 7 45 164 3 0 fail 
SANE 0.0 1 40 36.44 0.00 356 7 46 159 1 0 fail 
SANE 0.0 1 80 41.10 0.00 312 8 31 TI 5 1 fail 
SANE 0.0 1 160 48.91 0.00 238 gi 13 21 39 1 fail 
SANE 0.0 10 1 1.26 0.00 28 419 0 0 42 0 fail 
SANE 0.0 10 10 30.20 0.00 0 386 0 0 53 0 fail 
SANE 0.0 10 20 134.87 0.00 28 153 0 0 262 0 fail 
SANE 0.0 10 40 211.01 0.00 73 46 3 3 122 0 fail 
SANE 0.0 10 80 223.92 0.00 60 56 4 3 98 0 fail 
SANE 0.0 10 160 249.08 0.00 45 99 3 2 87 0 fail 
SANE 0.5 1 1 0.28 0.01 74 584 0 0 0 0 fail 
SANE 0.5 1 10 2.05 0.04 1 69 0 0 23 0 fail 
SANE 0.5 1 20 4.38 0.08 0 114 0 0 96 0 fail 
SANE 0.5 1 40 8.60 0.16 0 211 0 0 67 0 fail 
SANE 0.5 1 80 13.63 0.25 0 245 0 0 53 0 fail 
SANE 0.5 1 160 18.22 0.33 2 215 0 0 76 0 fail 
SANE 0.5 10 1 0.45 0.01 203 557 198 189 3 0 fail 
SANE 0.5 10 10 5.53 0.10 0 185 0 0 22 0 fail 
SANE 0.5 10 20 17.96 0.33 0 58 0 0 39 0 fail 
SANE 0.5 10 40 56.18 1.03 0 278 0 0 137 0 fail 
SANE 0.5 10 80 110.23 2.02 2 204 0 0 125 0 fail 
SANE 0.5 10 160 140.88 2.58 2 107 0 0 93 0 fail 
SANE 0.94 1 1 0.05 0.02 165 59 0 0 0 0 fail 
SANE 0.94 1 10 0.32 0.13 336 594 481 439 21 0 fail 
SANE 0.94 1 20 0.73 0.30 61 521 0 0 12 0 fail 
SANE 0.94 1 40 1.35 0.56 45 525 6 35 3 0 fail 
SANE 0.94 1 80 2.00 0.83 122 420 14 121 1 0 fail 
SANE 0.94 1 160 2.18 1.20 182 321 13 163 0 0 fail 
SANE 0.94 10 1 0.07 0.03 181 33 0 0 0 0 fail 
SANE 0.94 10 10 0.49 0.20 582 469 113 107 116 0 fail 
SANE 0.94 10 20 1.57 0.65 0 427 0 0 137 0 fail 
SANE 0.94 10 40 5.86 2.44 0 493 0 0 122 0 fail 
SANE 0.94 10 80 12.76 5.31 0 272 0 0 178 0 fail 
SANE 0.94 10 160 17.71 7.37. 1 144 0 0 205 0 fail 
MAD —0.94 1 1 0.13 1.74 118 23 2 41 109 0 fail 
MAD —0.94 1 10 0.20 2.60 112 116 0 3 8 0 fail 
MAD —0.94 1 20 0.24 3.15 94 199 0 0 0 0 fail 
MAD —0.94 1 40 0.31 3.98 92 256 0 7 0 0 fail 
MAD —0.94 1 80 0.41 5.31 134 261 4 51 1 0 fail 
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Table 3 
(Continued) 
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Note. The number of images passing each polarimetric constraint are given along with the number Mim simultaneously passing all of them. The accretion rate M. is 
in units of 10-7 Mo yr !, and Piet42 is the jet power in units of 10? erg ѕ!. Models that pass according to either the simultaneous or joint scoring method (boldface) 
and have Pjet42 2 1 are given a summary score of pass. 
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